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ABSTRACT 

The  objective  of  this  paper  is  to  describe  the  Naval  Surface 
Weapons  Center  (NSWC)  Program  for  developing  high  performance, 
simple,  rugged,  cost  effective  magnetoelastic  force  feedback 
sensors  for  robots  and  machine  tools.  Recent  advances  in  magneto¬ 
elastic  material  technology  have  paved  the  way  for  corresponding 
improvements  in  the  state  of  the  art  in  force  feedback  sensors 
for  robots  and  machine  tools.  Also,  NSWC  has  designed  magnetic 
circuits  which  are  easily  adapted  to  force  feedback  sensors.  In 
this  paper,  magnetoelastic  materials  are  described  along  with  the 
properties  that  make  them  potentially  such  outstanding  force  feed¬ 
back  sensors.  Following  this,  the  NSWC  Program  is  detailed  in¬ 
cluding  advances  in  materials  research,  in  simple,  low  cost 
electronic  and  magnetic  circuits,  and  designs  for  force  feedback 
sensor  modules.  The  results  are  in  the  public  domain. 

INTRODUCTION 

NSWC  is  developing  high  performance,  simple,  rugged,  cost  effec¬ 
tive  magnetoelastic/magnetostrictive  force  feedback  sensors  for 
machine  tools  and  robots.  (Note  figure  1.) 

The  Navy  is  facing  increasingly  severe  manpower  and  skills  short¬ 
ages.  Costs  of  equipment  and  complexities  of  equipment  mainte¬ 
nance  are  escalating.  To  combat  this,  CAD/CAM,  intelligent 
automation,  and  robotics  are  depended  upon  to  play  a  vital  role. 
Tactile/force  feedback  sensors  will  be  at  the  center  of  the  Navy 
effort . 

To  optimize  its  contribution,  NSWC  is  pulling  together  a  range  of 
disciplines,  technical  expertise,  and  experience  ranging  from 
materials  basic  research,  to  signal  processing  techniques,  to 
modular  robotic  sensor  designs. 

The  heart  of  this  coordinated  effort  is  the  materials  basic  re¬ 
search.  Since  the  1960's,  NSWC  has  been  performing  basic  research 
in  rare  earth  materials.  In  recent  months,  breakthroughs  have 
been  achieved  which  permit  these  materials  to  be  used  to  push  the 
state  of  the  art  for  a  range  of  practical,  high  performance, 
tactile/force  feedback  sensor  applications. 
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NSWC  has  also  made  major  technical  accomplishments  in  magnetometer 
circuitry  for  sensitive  magnetic  field  sensors.  This  expertise  is 
being  applied  toward  the  closely  related  problem  of  force  feedback 
sensor  circuitry  and  signal  processing. 

In  this  paper,  magnetoelastic/magnetostrictive  materials  are  de¬ 
scribed  along  with  the  properties  that  potentially  make  them  such 
outstanding  force  feedback  sensors.  Following  this,  the  magnetic 
circuit  signal  processing  techniques  will  be  described.  Finally, 
all  these  disciplines  will  be  integrated  to  form  tactile/force 
feedback  sensor  modules  for  grip  and  torque.  Designs,  circuits, 
and  research  results  are  in  the  public  domain. 

MATERIALS 

In  this  section,  the  magnetoelastic  effect  is  explained  including 
its  subsets,  the  "Villari"  effect,  and  the  magnetostrictive 
effect.  Next,  the  make  up  and  characteristics  of  the  NSWC  mate¬ 
rials  are  described  and  the  reasons  they  are  potentially  such 
outstanding  tactile/force  feedback  sensors.  Finally,  research 
test  results  are  given. 

A  material  is  magnetoelastic  if  there  is  a  relationship  between 

(1)  changes  in  the  internal  magnetic  moment  (and  hence  B  field), 

(2)  changes  in  the  physical  forces  applied  to  or  by  it,  and  (3) 
changes  in  its  physical  length.  A  magnetoelastic  material  ex¬ 
hibits  the  "Villari"  effect  when  it  shows  a  change  in  its  B 
field  as  a  result  of  its  being  subjected  to  external  forces  of 
tension  or  compression.  A  magnetoelastic  material  exhibits  the 
magnetostrictive  effect  when  it  shows  a  change  in  either  (1)  its 
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physical  length  due  to  changes  which  have  been  induced  in  its 
internal  B  field,  or  (2)  the  inverse-changes  in  its  internal 
B  field  which  have  been  induced  by  changes  in  its  physical 
length . 

In  this  paper  we  are  basically  concerned  with  two  families  of 
magnetoelastic  materials;  amorphous  ribbons  -  Fe70CoioB2o  (with 
small  amounts  of  silicon  and  cobalt)  and  Tb . 2 7Dy . 7 3Fe2  rods. 

Both  materials  exhibit  "Villari"  and  magnetostrictive  effects.* 

The  amorphous  ribbons  are  superior  in  "Villari"  effect  sensing 
applications  due  to  stretching  forces  and  the  Tb . 2 7Dy . 7 3Fe2  rods 
are  superior  in  magnetostriction  (but  they  can  also  perform 
excellent  "Villari"  effect  sensing  in  response  to  compressive 
forces ) . 

This  materials  development  was  begun  when  NSWC  scientists  reasoned 
that  while  rare  earth  materials  possess  extraordinary  magnetic 
properties  at  below  room  temperatures,  it  might  be  possible  to 
develop  alloys  and  amorphous  solutions  combining  the  rare  earths 
with  the  more  classical  ferrous  materials  to  provide  practical 
new  materials  with  extraordinary  performance  properties.  These 
materials  can  be  used  at  normal  room  temperature  and  above.  After 
nearly  20  years  of  research  into  this  matter,  the  materials  are 
now  reaching  the  point  where  they  are  ready  for  industrial 
applications . 

Let  us  now  discuss  how  these  materials  act  as  tactile/force  feed¬ 
back  sensors.  As  shown  in  figure  2,  the  materials  are  first 
treated  so  that  each  internal  magnetic  domain  lines  up  with  its 
net  magnetic  moment  perpendicular  to  the  long  axis  of  the  material. 
At  the  same  time,  the  net  magnetic  moment  of  each  domain  is  pointed 
in  a  direction  opposite  to  that  of  its  neighbor's.  This  leaves  the 
material  with  a  total  net  magnetic  moment  of  zero;  thus  reducing 
spurious  effects  and  increasing  material  sensitivity.  If  the 
magnet^elastic  material  is  now  biased  by  an  external  magnetic 
field  H  (either  by  permanent  magnetic  material  or  current)  the 
magnetic  moments  rotate  toward  the  direction  of  the  field.  For 
maximum  linear  dynamic  range,  the  bias  H  field  is  set  such  that 
the  direction  of  the  net  magnetic  moment  of  each  magnetic  domain 
is  at  approximately  45°  to  the  axis  of  the  material. 


*Magnetoelastic  materials  research  has  been  in  process  at  NSWC 
since  the  1960's.  The  material  "terfenol"  (Tb . 2 7Dy  .  7 3Fe2 )  was 
discovered  by  Dr.  A.  E.  Clark  at  NSWC.  Dr.  H.  Savage  and  staff 
have  pursued  the  work  in  perfecting  the  material  (along  with 
several  applications  patents).  Research  in  another  material 
(amorphous  ribbons)  was  begun  at  NSWC  by  Dr.  M. Mitchell  and  has 
since  been  continued  by  Dr.  Kabacoff  and  staff. 
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•  IF  WE  INTEGRATE  WITH  RESPECT  TO  TIME 
WE  HAVE  A  MEASURE  FOR  FORCE 


FIGURE  2  HOW  THE  MATERIALS  ACT  AS  FORCE  SENSORS 


Magnetoel asticity 

As  the  ribbon  is  stretched  by  an  external  force,  the  net  magnetic 
moment  of  each  domain  rotates  toward  the  ribbon  axis  at  an  angle 
which  is  linearly  proportional  to  the  exciting  force.  The  addi¬ 
tive  effect  of  all  these  magnetic  moments  rotating  toward  the 
ribbon  axis  is  that  the  ribbon  total  B  field  rotates  toward  its 
longitudinal  axis,  linearly  proportional  to  the  exciting  force. 

By  the  same  token,  as  a  terfenol  rod  is  compressed  by  an  external 
force,  its  total  B  field  rotates  away  from  its  longitudinal  axis, 
linearly  proportional  to  the  exciting  force.  In  both  cases,  this 
is  a  "Villari"  effect.  This  linear  relationship  can  be  described 
by  the  following  equations: 


where  |Ab|  is  the  total  net  change  in  the  Material  B  field,  d  is 

AT?  I 

a  constant  and  |-j'  is  the  applied  force  per  unit  area. 
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where  y  is  the  Young’s  Modulus  of  the  material  in  psi  and  |^|-| 
is  the  material  change  in  length  per  unit  length. 

=  d [ AH |  where  | AH |  is  an  externally  applied  magnetic  field. 

(One  should  recall  that  §  =  p  3  where  y  is  the  material 
permeability  (uQii  ). 

Referring  to  figure  2,  if  we  wrap  either  an  amorphous  ribbon  or  a 
terfenol  rod  with  coils  of  insulated  wire  and  send  current  through 
the  wire,  a  $  field  will  be  superimposed  on  the  B  field  already  in 
the  material.  This  will  cause  the  amorphous  ribbon  and  terfenol 
rods  to  expand  or  contract  depending  on  the  product  of  the  cur¬ 
rent  and  the  number  of  coils.  This  magnetostrictive  property  is 
most  pronounced  for  terfenol  rods  (see  figure  3). 


CLOSED  PATH  OF 
INTEGRATION 


NOT  TO  SCALE 

FIGURE  3  H&B  FIELDS  IN  THE  MAGNETOSTRICTIVE 
ROD 


Again,  we  can  use  the  |AB|  =  d | AP/A |  relationship.  But  A$  =  yAft 
where  y  =  p  Pr  =  4  x  10“7(8)  since  y  terfenol  =  8.** 


The  question  then  becomes  how  large  can  we  make  AB?  Let  us  illus¬ 
trate  by  picking  a  few  numbers.  Using  a  current  of  1  amp,  a  rod 
length  of  2  inches  (5.08  cm),  .25  in.  diameter  (.635  cm)  and  with 
2000  turns  of  wire  around  it,  we  will  get: 

8  •  dt  =  I  =  2000  (1) 


(where  length  is  in  meters).  And  referring  back  to  figure  3,  this 
simplifies  to  A b  (.0508)=  2000  (1  amp) 


Thus  5  =  .126  webers/m2  j 

| b I  =  d | AF/A | ;  d  =  .25  x  lO"7^ 
1  lb  =  4.44  newtons.  / 

So  $  =  35.9  lbs  or  159  newtons 


webers/newton  for  terfenol  and 


Again  back  to  our  equation 

|  AS  |  =  d|^-|  wheri  d  =  ^  ^  x  102  webers/m2 

o  r 

Ab  =  8.6  x  102,  f^ln 
But  AB=  .126  webers/m2 

Thus  ’Kl  =  2.93  x  10~4  in  ( 7  /  4  4  x  10-4  cm) 

Magnetostriction 

Of  course  one  can  see  that ; increasing  the  product  of  the  current 
and  windings  increases  bot.'  the  force  and  the  magnetostrictive 
travel.  (Magnetostrictior  ;  of  .001  in.  per  inch  is  easily 
achievable  for  terfenol.)  ■ 


Force  Sensing 

There  are  essentially  tw/  methods  for  force  sensing  using  magneto¬ 
elastic  materials,  one  ?/ technique  which  measures  the  total  force 
acting  on  the  sensor  an</  one  which  measures  the  rate  of  change  of 
the  force  acting  on  the  'sensor. 


i 


**Values  and  equations  .  rom  Dr.  Howard  Savage’s  notes  4/13/81. 
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First  we  will  discuss  the  technique  which  measures  the  total  force 
on  the  sensor.  In  this  technique  an  external  oscillator  source  is 
used.  It  should  be  noted  that  the  amorphous  ribbon  with  its  y  of 
the  order  of  magnitude  of  20,000  (vs  8  for  the  terfenol  rods)ris 
the  better  of  the  two  for  the  total  force  sensing  technique.  So, 
for  this  method  we  shall  use  the  amorphous  ribbon  in  our  example. 
Let  us  start  by  coiling  an  insulated  wire  around  an  amorphous 
ribbon  and  connecting  this  wire  to  an  AC  source.  This  will 
provide  the  drive  circuit.  (See  figure  4.)  We  can  wrap  another 
set  of  insulated  wire  coils  over  the  coils  of  the  driving  circuit 
and  these  provide  the  signal.  The  driving  circuit  is  constructed 
so  that  it  drives  the  amorphous  ribbon  back  and  forth  into 
saturation.  Since  the  waveform  drops  back  rapidly  to  zero,  the 
output  signal,  which  is  proportional  to  d$,  is  maximum.  But  we 

dt 

know  that  the  characteristics  of  the  amorphous  ribbon  are  such 
that  its  ur  changes  with  applied  force.  This  means  that  the  sat¬ 
uration  B  field  changes  and  we  have  a  difference  between  the  B 
field  with  and  without  applied  force  (Ae).  This  Ab  is  some  func¬ 
tion  of  the  applied  force  or  pressure  which  we  can  call  fi  (P). 
This  Ab  translates  into  a  difference  in  the  voltage  output  V  which 
in  turn  relates  back  to  the  applied  force  or  pressure  by  some 
function  F2(P).  We  call  this  the  external  oscillator  technique. 

( See  figure  5 . ) 


At  =  I2|P) 
AV=  fi(P) 


FIGURE  4  SIGNAL  PROCESSING  SIN  MAGNETOELASTIC 
RIBBONS 
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FIGURE  5  THE  SEPARATELY  EXCITED  OSCILLATOR 


Next  we  will  discuss  the  technique  which  measures  the  derivative 

/-5  Li  I 


of  the  total  force  (or  pressure)  with  respect  to  time  (g^r)  • 
technique  depends  on  the  rate  of  change  of  flux  <j> 

(<j)  =  £  g  •  ds  (s  =  area  of  ribbon  or  rod  cross  section). 

governing  equation  is: 


This 


=  n  ^ 
dt 


n  =  number  of  coils 
V  =  voltage  output 


Looking  at  figure  2,  we  can  see  that  we  have  essentially  a  passive 
sensor  with  one  set  of  coils  (no  driving  circuit).  It  will  oper¬ 
ate  at  very  low  frequencies  (to  DC) ,  those  of  the  robotic  tactile 
force  itself.  At  these  low  frequencies,  and  with  the  small  cross 
sectional  areas  of  the  ribbons  and  rods  (typically  .006  in2),  there 
is  virtually  no  skin  effect  or  magnetic  hysteresis  and  the 
equation  simplifies  to: 

V  =  n  f  A3  •  ds  (AB  =  d^) 


(S  =  A) 


„  .  dF 

V  =  nd2  dt 


Thus  we  can  see  that  the  voltage  is  proportional  to  the  term 
(or  |gf*|)  and  the  integral  of  J ^  V  dt  yields  the  total  force. 

The  high  efficiency  of  the  rods  and  ribbons  (mechanical  work  on 
the  material  compared  to  the  change  in  magnetic  field  Ab)  makes 
the  derivative  sensor  practical. 
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Generally  the  force  derivative  sensor  provides  a  relatively  low 
voltage,  high  power  output.  And,  comparing  figures  2  and  3,  one 
can  see  that  when  a  terfenol  rod  is  used  there  is  the  option  of  it 
acting  either  as  a  sensor  or  actuator.  On  the  other  hand,  the 
external  oscillator  technique  yields  higher  output  voltages  and 
the  sampling  time  is  not  critical  (since  total  force  is  measured). 
Clearly,  the  preferable  technique  depends  on  the  application. 

Let  us  now  discuss  why  these  materials  are  potentially  such  good 
force  feedback  and  tactile  sensors. 

Outstanding  Sensitivity 

Both  the  ribbons  and  the  terfenol  rods  have  demonstrated  outstand¬ 
ing  sensitivity.  A  0.8  ratio  for  magnetomechanical  coupling  has 
already  been  produced  at  NSWC  for  both  the  rods  and  amorphous  rib¬ 
bon  and  the  value  is  expected  to  climb  to  0.9.  This  0.8  value  is 
the  ratio  of  the  output  voltage  to  the  voltage  equivalent  of  the 
force  applied  to  the  material.  A  40  mV/V  ratio  was  previously 
considered  excellent. 

Outstanding  Dynamic  Range 

A  conservative  estimate  of  the  ribbon's  linear  region  for  dynamic 
range  is  +1,000  psi.  It  should  also  be  noted  that  experiments  at 
NSWC  (for  underwater  pressure  sensors)  have  indicated  that  it  is 
possible  to  resolve  .004  psi  in  a  background  depth  of  266  psi. 

This  means  that  theoretically  we  can  expect  10  log  (266/. 004)  = 
48.2  dB  volts.  30  dB  power  or  15  dB  volts  is  normally  considered 
to  be  excellent.  The  terfenol  rods  have  a  dynamic  range  2  times 
that  of  the  ribbons  and  are  prestressed  to  2,000  psi. 

No  Observable  Mechanical  Hysteresis 

The  amorphous  ribbons  and  terfenol  rods  are,  for  all  intents  and 
purposes,  completely  recoverable.  This  recoverability  extends  to 
well  beyond  2  or  3  times  the  +1,000  psi  linear  dynamic  range  of 
each. 

Outstanding  Linearity 

Graphs  for  both  experiments  conducted  by  Japanese  investigators, 

K.  Mohri  and  E.  Sudoh,  showed  outstanding  linear  response  and 
this  linearity  has  been  duplicated  at  NSWC  by  investigator 
J.  P.  Scarzello . *** 

Simple  Electronics 


The  circuitry  for  the  oscillator  drive  method  is  shown  in  figure  5 
As  can  be  seen  from  the  figure,  only  one  small  circuit  board  is 
needed . 


***P.4  NSWC  memo  dtd  23  Nov  19 7 9 >  "Development  of  a  pressure  trans 
ducer  using  amorphous  magnetic  materials,"  by  J.  F.  Scarzello. 
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Tough,  Inexpensive,  Corrosion  Resistant  and  Radiation  Hardened 
Materials 


The  amorphous  ribbons  are  tough  and  can  stand  up  to  20,000  psi 
(13,764  newtons/cm2).  The  rods,  when  prestressed  to  2,000  psi 
(1,376  newtons/cm2)  in  a  metal  can,  are  also  tough.  Both  materials 
can  be  mass  produced  inexpensively.  Both  are  extremely  resistant 
to  corrosion  and  radiation.  But  there  are  aspects  of  the  material 
which  require  additional  investigation  and/or  are  a  bit  troublesome. 

Stray  Fields  and  Voltage  Offsets 

This  is  a  problem,  particularly  for  the  amorphous  ribbons  because 
they  have  such  a  high  y  (magnetic  permeability)  and  because  they 
are  driven  by  the  "external  oscillator"  technique.  NSWC  is 
oursuing  three  techniques  to  deal  with  this  problem:  (1)  putting 
magnetic  shielding  material  around  the  ribbons,  (2)  using  2  sets 
of  symmetrical  windings  so  that  the  errors  are  self  cancelling, 
and  (3)  using  filtering  and  digital  signal  processing  techniques. 

Noise  in  the  Ribbon  (when  operated  as  part  of  the  external  drive 
oscillator  system) 

This  noise  tends  to  restrict  the  available  dynamic  range  of  the . 
system.  As  a  result,  the  materials  experts  at  NSWC  are  continuing 
to  work  on  improving  the  material  and  considerable  progress  has 
been  made.  It  has  also  been  noted  that  the  ribbons  tend  to 
resonate  magnetostrictively ,  so  selecting  the  proper  drive 
frequency  can  lower  noise  considerably.  15  KHz  has  yielded  good 
results . 

Brittleness  in  the  Terfenol  Rods 


Very  recent  experiments  at  NSWC  have  resulted  in  terfenol  rods 
which  can  stand  2,000  psi  (1,376  newtons/cm2)  in  tension  (in 
addition  to  high  compressive  capabilities).  This  represents 
still  another  improvement  in  the  breed. 

Low  ur  in  the  Terfenol  Rods 

Low  ur  and  high  magnetostrictive  capabilities  go  hand  in  hand  so 
only  limited  progress  can  be  made  in  this  area.  This  means  that 
the"  terfenol  rods  must  act  as  derivative  sensors  and  cannot  be 
driven  by  the  separately  excited  oscillator  technique. 

Temperature  Stability 

This  has  not  been  fully  tested  though  stability  to  100°C  is 
anticipated . 
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SIGNAL  PROCESSING 


In  this  section,  signal  processing  techniques  will  be  outlined 
including  both  the  oscillator  drive  technique  and  the  derivative 
method.  Some  of  the  expected  signal  processing  problems  and 
their  usual  methods  of  solution  are  addressed. 

The  oscillator  drive  method  is  the  main  signal  processing  tech¬ 
nique  to  be  used  in  the  grip,  torque,  and  grip/torque  modules. 

As  mentioned  above,  this  method  measures  the  total  force  or  torque 
as  the  case  may  be  and  it  is  simple,  sensitive,  and  accurate.  In 
addition,  this  technique  is  well  understood  at  NSWC  since  it  was 
developed  for  magnetometers  in  underwater  mine  and  space  applica¬ 
tions,  and  it  has  resulted  in  many  major  publications  and  patents 
dating  back  to  1970.  (For  example,  a  magnetometer  designed  jointly 
by  NSWC  and  NASA  Ames  Research  Center  was  a  part  of  the  equipment 
deployed  on  the  lunar  surface  in  the  Apollo  16  Mission.) 

The  Relationship  of  Magnetometers  to  Robotic  Force  Feedback  Sensors 

In  a  magnetometer,  the  magnetically  active  element  (in  our  case  an 
amorphous  ribbon)  has  its  B  field  changed  by  an  intruding  magnetic 
field.  This  B  in  turn  changes  the  magnetic  reluctance  of  the 
amorphous  material.  If  the  material  is  being  excited  by  an  oscil¬ 
lator  drive,  as  an  intruding  magnetic  field  adds  to  or  subtracts 
from  the  magnetic  field  in  the  magnetometer  core,  output  voltage 
will  be  affected.  But  as  we  noticed  earlier,  the  same  effect  will 
occur  if  the  ribbon  is  physically  stressed.  Thus  the  oscillator 
drive  signal  processing  technique  is  equally  applicable  to  force 
feedback  sensors  and  magnetometers. 

A  Brief  Description  of  how  the  Oscillator  Drive  Processing 
Technique  Works 


A  key  factor  is  the  current/field  wave  shape  in  the  driving  cir¬ 
cuit  coils.  This  wave  shape  is  such  that  the  current/field  builds 
to  a  peak  and  then  rapidly  drops  to  zero.  The  resulting  large 

^  can  be  used  to  provide  a  voltage  output  which  is  a  function  of 

the  force  on  the  ribbon  (see  figures  4  and  5). 

As  the  sensor  ribbon  is  stretched  by  a  robotic  force  or  torque, 
the  slope  of  the  B-^  curve  becomes  steeper  and  the  output  voltage 
greater.  Also,  inductance  is  proportional  to  the  slope  of  this 
curve,  so  a  circuit  which  will  produce  a  signal  proportional  to 
inductance  will  produce  a  measurable  signal. 

The  central  component  of  the  circuit  is  an  operational  amplifier 
(op-amp).  As  a  result  of  high  amplifier  gain.  A,  and  high  input 
impedance,  it  can  be  shown  that  the  output  voltage  is  related  to 
the  input  voltage  by  the  ratio  of  feedback  and  input  impedance. 
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Suppose  the  input  voltage  is  sinusoidally  varying: 

V  =  V  sin(wt) 
a 

Its  magnitude  is  given  by 
|Val  =  V 

If  the  input  component  is  a  resistor,  and  the  feedback  component 
an  inductor:  then 

Za  ■  R  Zb  -  J“L 


and 


M  -  <^1  t 


Thus  a  signal  is  produced  which  is  proportional  to  inductance. 

If  both  ribbons  are  pre-tensioned ,  an  additional  force  will 
increase  inductance  in  one  coil,  and  de c rease  inductance  in  the 
other,  both  by  L  (figure  5). 

The  voltage  magnitudes  at  points  and  V2  are  given  by 

Tr  _  Vo)  (L  -AL)  „  _  Vm  (L  +AL) 

V1 - R  V2  R 

The  orientation  of  the  diodes  cause  the  dc  levels 

V3  -  |VX|  v„  =  -|V2| 

The  final  op-amp  adds  and  V^. 

V.  =  -(1^)  (V,  +  V,) 


O 


Ra' 


Combining  the  above 


V  —  ( ^b  wt 

Vo  ”  (— R"R^)AL 
a 
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This  is  a  signal  proportional  to  the  difference  in  inductance 
caused  by  an  imbalance  in  ribbon  tension  (sensor  ribbon  compared 
to  a  ref.  ribbon).  This  simplified  circuit  is  expanded  to  include 
a  zero  adjustment  and  ideal  diodes. t  Specific  component  values 
are  shown  in  figure  5 . 

The  V0  from  the  circuit  shown  in  figure  5  can  now  be  further 
processed.  We  can  pick  up  the  RMS  voltage  output  (using  perhaps 
a  simple  square  law  detector).  In  this  process,  the  equation 
VQ  =  4.44  N2f  BA  x  10~8  volts 

N  =  number  of  windings 
f  =  frequency 
A  =  area 

is  a  good  estimate  of  the  RMS  value  we  would  expect.  Since  a 
force  on  a  robotic  sensor  introduces  a  change  in  B  using,  once 

again,  the  equation  |ab|  =  d|jp| ,  we  will  get  a  change  in  the 

total  voltage  output  of  approximately: 

V  =  4.44  N~f  B  A  x  10“8  volts 
rms  d 

since 


AV  xlO-8 

I  ao  |  =  rms _ 

* AF|  4.44  N2fd 

Thus  simply  measuring  the  RMS  voltage  output  will  provide  us  with 
the  information  we  need  to  solve  the  force,  or  torque,  directly. 

The  derivative  technique  can  also  be  used  for  signal  processing  in 
any  of  the  magnetoelastic/magnetostrictive  robotic  force  feedback 
sensors.  As  mentioned  above,  in  the  derivative  technique  the 

equation  V  =  n  ^  applies  or  V  =  n  ^  j  dAF  | 

V  =  n  d|^r|.  So,  again,  we  are  dealing  with  the  rate  of  force 
application . 


tFrom  "Progress  Report  Magnetoelastic  Sensor  Development,"  6/15/81 
thru  8/15/81  Prof.  R.  DeMoyer  and  Prof.  E.  Mitchell,  U.S.  Naval 
Academy,  9/15/81,  pp.  11-13 
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ROBOTIC  TORQUE,  GRIP,  GRIP/TORQUE  MODULES 

In  this  section,  preliminary  designs  for  robotic  torque  and  grip 
modules  will  be  described.  In  addition,  since  the  sensors  are 
designed  to  be  modular,  a  grip/torque  sensor  will  be  described 
which  is  the  result  of  cascading  a  grip  and  a  torque  sensor  module 
together.  The  modules  are  for  a  large  industrial  robot  that  can 
lift  250  lbs  (114  Kg) . 

The  Torque  Module  (figures  4,  5  and  6) 

The  design  goals  for  the  torque  module  are  125  in-lb  (144.1  Kgcm) 
of  torque  in  the  linear  region  in  a  package  of  nominally  2.5  in. 
diameter  (6.35  cm)  and  3/4  in.  (1.9  cm)  height.  Again,  as  shown 
in  figure  1,  this  module  is  designed  to  go  into  the  wrist  of  the 
robot  gripper.  The  125  in-lb  torque  parameter  is  picked  because 
it  is  assumed  that  a  worst  case  condition  is  one  in  which  the 
robot  picks  up  a  225  lb  object  with  a  .5  in.  deviation  from  the 
object’s  center  of  gravity.  With  0.030  in.  (.08  cm)  typical 
robot  movement  accuracy,  this  represents  an  extreme  case. 

Design  Calculations 

Figure  6  shows  the  design  philosophy.  When  the  flat  bar  spring 
bends,  the  ribbon  is  either  stretched  or  compressed  (depending  on 
direction  of  torque).  This  is  a  linear  relationship  with  the  force 
on  the  end  of  the  bar.  From  our  design  requirements  we  know  that 
we  need  a  maximum  torque  of  125  in-lbs.  This  must  be  balanced  by 
the  bending  moment  in  the  aluminum  bars  of  the  torque  sensor. 

/•X 

1^1  (125  in-lbs)  =  ( 4 ) ( 2 )  /  x°dF 

Jo 

4  =  number  of  bars 
2  =  tension  and  compression 
xdF  =  Atorque 

Fx  C 

For  a  rectangular  beam  J  wdx  =  dF 

|^-|  =  psi  at  outer  fiber  of  bar 
x  =  distance  from  bar  center  to  outer  fiber 


But 


w  =  bar  width 


Y  =  10.3  x  106  psi 
m 

for  aluminum 
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SENSOR  HOUSING 
REACTIVE 
FORCES 


FORCE 

FIGURE  6  TORQUE  SENSOR 


and  for  amorphous  ribbons  =  25  x  10  6in/in 
Lets  add  on  25  x  10~6in/in  for  adhesive  shear. 


If  we  make  it  3/4  in.  wide,  it 
that  if  our  electronics  has  30 
.125  in-lb  to  125  in-lb  in  the 
linear  region.  Of  course,  we 


must  be  .4  in .  thick.  This  means 
dB  dynamic  range,  we  can  pick  up 
linear  to  250  in-lb  using  the  non¬ 
will  expect  much  better  sensitivity. 
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We  know  V. 


=  4.44  N9f  BA  x  10“8  volts 
rms  d 

B  =  d  P/A 


(CGS  units) 


for  Mks  units  use 


d 


11.1 


Plb  tt 
2 

Acm 


with  f  =  15  KHz;  V  _  required  =  . lmV 

rms 


we  get  the  design  relationship 

.0135  =  P  lb  N,  or  .00614  =  P,  N0 

^  kg  2 

We  will  use  a  ribbon  .2  in.  wide  by  .002  in.  thick  with  0  to 
+1,000  psi  linear  dynamic  range.  So,  at  its  maximum  stress,  it 
will  experience  (1,000)  ,002(.2)  =  .4  lb  tension  or  compression. 
This  corresponds  to  .125  in-lbs  torque  on  the  beams.  To  get  down 
to  .025  in-lbs  (28.9  gmcm  or  37  dBV  dynamic  range),  we  need 
.00008  lb  sensitivity  in  the  ribbon: 

.00008  N2  =  .0135 


N2  =  168;  use  N2  =  170  turns 

It  will  be  a  challenge  to  get  the  electronics  to  handle  such 
sensitivity . 

At  this  point  we  will  examine  the  results  from  a  test  prototype 
torque  sensor.  This  prototype  was  built  oversize  so  as  to  facili¬ 
tate  testing  and  modification.  It  is  essentially  patterned  on  the 
concept  shown  in  figure  6. 

Results  (see  figure  7) 

At  10  KHz,  inductance  increases  with  tension  up  to  a  point  after 
which  it  remains  constant.  This  is  true  for  both  annealed  and 
unannealed  ribbon. 

Slight  deviations  from  linearity  at  the  curve  extremes  are 
probably  due  to  the  adhesive  yield  in  shear. 

Using  annealed  ribbon,  the  change  in  inductance  effect  with  tension 
is  reversed  from  the  unannealed  version. 


ttMeeting  3/23/80,  Dr.  Howard  Savage,  J.  M.  Vranish. 
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The  annealed  ribbon  can  take  high  signal  levels  without  saturation, 
and  the  output  is  a  reasonably  linear  function  of  applied  load  in 
both  tension  and  compression. 

Conclusions 


The  annealed  ribbon  excited  at  a  frequency  near  its  resonance 
exhibits  desirable  properties:  near  linearity,  large  dynamic  range, 
and  high  level  output.  The  phase  shift  between  the  annealed  and 
unannealed  is  due  to  the  fact  that  annealed  and  unannealed  ribbons 
resonate  magnetostrictively  at  different  frequencies. 

Future  Work 


A  magnetic  return  path  should  be  installed  to  determine  if  the 
effects  of  external  iron  proximity  can  be  totally  eliminated. 

Many  trade-offs  include: 

Signal  levels 
Frequency 

Number  of  turns  on  coils 

All  components  must  be  miniaturized  in  order  to  construct  a 
practical  torque/force  sensor. tit 

We  can  now  show  the  design  approach  that  will  be  taken  on  the  grip 
sensor.  This  is  the  same  as  the  torque  sensor  except  the  force 
comes  straight  down  and  the  ribbons  are  on  the  top  and  bottom  of 
the  bars  rather  than  the  sides  as  in  the  torque  sensor. 

The  design  goals  for  the  grip  sensor  are  250  lbs  (182  kg)  maximum 
grip  in  a  small  package  2.5  in.  (6.35  cm)  diameter  and  1  in. 

(2.5^  cm)  deep. 

Grip  Sensor  Design  Calculations 

The  grip  sensor  will  be  designed  similar  to  the  torque  sensor. 

250  lbs  (113.5  kg)  maximum  design  grip  will  be  distributed  among 
the  4  aluminum  bars  of  the  grip  sensor. 


So:  62. 5#-  =  8  Jx  -  -  wdx  as  before.  Z  =  bar  length 

*0  o 

and  this  reduces  to  62. 5&  =  4118  wx  ^ 


tt+tiprogress  Report  Magnetoelastic  Sensor  Development,"  6/16/81 
through  8/15/81,  pp.  28-31. 
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using  i  =  1  in.  and  picking  w  =  3/8  in.,  we  get  a  bar  thickness 
of  .5  in.  (1.02  cm).  Our  electronics  relationships  remain  the 
same  as  for  the  torque  sensor. 

.0135  =  Flb  N2 

If  we  use  .008  lb  in  our  ribbon  as  the  minimum  force  to  pick  up 
at  . lmV  rms  signal  output  voltage,  we  will  need  170  coils  as 
before . 

Future  Directions 

During  the  first  year,  NSWC  will  continue  to  develop  and  refine 
the  sensor  modules.  This  will  include  optimizing  the  amorphous 
sensor  material  for  robotic  applications. 

In  the  second  year,  the  sensor  modules  will  be  interfaced  with  the 
NSWC  robot  to  iron  out  the  control  theory  and  vision  coordination 
questions.  Also  during  this  year,  a  gripper  will  be  designed  for 
the  removal  of  high  tension  fasteners  from  Navy  missiles. 

In  the  third  year,  a  gripper  will  be  built  and  the  sensor  modules 
interfaced . 

In  the  fourth  year,  the  gripper  will  be  interfaced  to  a  heavy  duty 
industrial  robot  and  the  combined  system  tested  in  a  Naval  Weapons 
Station. 

SUMMARY  AND  CONCLUSIONS 

In  summary,  we  have  described  the  NSWC  Program  for  developing  high 
performance,  simple,  rugged,  cost  effective  magnetoelastic/ 
magnetostrictive  force  feedback  sensors  for  machine  tools  and 
robots  in  CAD/CAM  operations.  We  have  shown  that  the  NSWC  Program 
is  a  comprehensive  program  including  basic  materials  research, 
signal  processing,  and  robot  sensor  modules.  We  have  outlined  the 
materials  characteristics,  the  signal  processing  techniques,  and 
the  robotic  sensor  designs.  The  basic  research  has  been  success¬ 
fully  completed,  and  practical  force  feedback  sensors  for  robots 
are  being  constructed  and  debugged.  The  NSWC  Program  will  signifi¬ 
cantly  advance  the  state  of  the  art  in  force  feedback  sensors. 
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ANNEALED  RIBBON 


TORQUE  (INCH-LB) 
FIGURE  7  LINEARITY  RESULTS 
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MODERN  CONTROL  TECHNIQUES  APPLICABLE  TO 
THE  SPACE  SHUTTLE  MAIN  ENGINE 


T.  C.  EVATT 


ROCKWELL  INTERNATIONAL  CORPORATION 
ROCKETDYNE  DIVISION 
SYSTEMS  ANALYSIS 
CANOGA  PARK,  CALIFORNIA  91304 


INTRODUCTION 


Historically,  development  of  liquid  propulsion  rocket  control 
systems  has  utilized  classical  control  techniques  to  regulate 
thrust  and  mixture  ratio.  The  approach  taken  here  is  to  model 
the  Space  Shuttle  Main  Engine  as  a  linear-multivariable  system 
whose  parameters  vary  with  engine  operating  environment.  Only 
inputs  and  noise  corrupted  outputs  realizable  at  an  actual 
engine  are  considered.  The  engine  control  objective  centered 
on  controlling  preburner  temperature,  mixture  ratio  and  thrust 
(combustion  pressure).  An  eighth  order  linear  model  with  two 
actuators  (valves)  and  their  associated  dynamics  is  derived  to 
apply  optimal  linear  quadratic  design  methodologies  to  control 
fuel  turbine  inlet  temperature  gradients.  The  design  metho¬ 
dology  selected  to  investigate  and  develop  a  feedback  con¬ 
troller  for  these  temperature  gradients  is  the  Linear  Quadratic 
Gaussian  (LQG)  design  philosophy.1  This  method  is  a  syste¬ 
matic  method  of  regulator  control  using  a  Kalman-Bucy  filter 
(state  estimator)  to  determine  plant  states  from  measured  para¬ 
meters.  A  subset  of  the  LQG  design  methodology,  the  Linear 
Quadratic  Regulator  (LQR)  will  be  derived  for  an  SSME  deter¬ 
ministic  model  at  one  design  point  and  graphical  results 
presented. 

The  use  of  this  modern  control  methodology  represents  an 
advance  in  the  design  philosophies  used  in  rocket  engine  con¬ 
trol  systems.  At  the  time  of  the  SSME  control  system  design  in 
the  early  1970's,  attention  was  given  to  both  classical  and 
modern  control  system  design  methodologies.  The  state  of  the 
art,  however,  was  not  mature  enough  to  use  the  modern  control 
systems  design  approach  for  the  SSME.  As  a  result,  an  advanced 
application  of  the  classical  methodology  was  selected  for  the 
control  system  design. 

With  the  foreseeable  increasing  performance  demands  on  the 
SSME,  LQG  control  methodology  that  relies  on  linear  quadratic 
synthesis  of  regulators  at  different  operating  points  is  being 
applied  to  design  an  improved  control  system.  Many  papers  have 
been  written  on  engine  control  using  Linear  Quadratic  Design 
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techniques.  Hoff  and  Hall2’3  have  discussed  a  control  design 
methodology  for  turbine  engine  control  synthesis  using  regu¬ 
lators  at  a  number  of  power  levels  at  sea  level,  static  condi¬ 
tions.  Taiwo1*  discussed  a  method  of  turbine  controller 
design  using  Zakain's  method  of  inequalities.  Merrill5  has 
used  sampled  data  regulators  to  design  multivariable  control 
laws  for  jet  propulsion.  Michael  and  Farrar5  discussed  a 
practical  and  systematic  linear  quadratic  synthesis  procedure. 
The  regulator  design  techniques,  however,  rely  on  full-state 
availability. 

This  paper  is  divided  into  four  sections.  The  first  section 
introduces  the  SSME  rocket  engine  and  the  techniques  used  to 
simulate  the  dynamics  of  the  SSME  including  linear  model  deri¬ 
vation.  A  discussion  of  actuator  dynamics  and  their  effect  on 
closed-loop  response  is  also  presented.  In  the  next  section, 
the  LQG  technique  is  briefly  discussed  with  a  more  detailed 
description  of  the  Linear  Quadratic  Regulator.  A  procedure  is 
discussed  for  choosing  the  LQR  weighting  matrices.  Conclusions 
and  recommendations  follow. 


SSME  TURBOPUMP  ROCKET  ENGINE 


INTRODUCTION 

With  the  advent  of  the  space  shuttle  concept  in  the  late  60’s, 
it  became  apparent  that  a  reusable  and  reliable  engine  design 
was  needed.  The  high  performance  requirements  of  the  Space 
Shuttle  Orbiter  demanded  a  new  state  of  the  art  in  liquid  pro¬ 
pellant  rocket  engine  design.  The  Space  Shuttle  Main  Engine 
(SSME)  design  uses  a  staged  combustion  power  cycle  (Figure  1) 
coupled  with  high  combustion  chamber  pressures.  In  the  SSME 
staged  combustion  power  cycle,  the  propellants  are  partially 
burned  at  low  mixture  ratio,  very  high  pressure,  and  relatively 
low  temperatures  in  the  preburners  to  produce  hydrogen-rich  gas 
to  power  the  turbopumps.  The  hydrogen-rich  gas  steam  is  then 
routed  to  the  main  injector  where  it  is  injected,  along  with 
additional  oxidizer  and  fuel,  into  the  main  combustion  chamber 
at  a  higher  mixture  ratio  and  pressure.  Liquid  hydrogen  is 
used  to  cool  all  combustion  devices  directly  exposed  to  contact 
with  high-temperature  combustion  products.  An  electronic 
engine  controller  automatically  performs  checkout,  start,  main- 
stage  and  engine  shutdown  functions. 

The  SSME  was  developed  especially  for  the  Space  Shuttle  Orbiter 
Vehicle,  which  uses  three  engines  for  launch.  The  SSME  is  a 
reusable,  high  performance,  liquid  propellant  rocket  engine 
with  variable  thrust.  The  engine  is  ignited,  on  the  ground  at 
launch  and  operated  in  parallel  with  the  solid  rocket  boosters 
during  the  initial  ascent  phase  and  continues  to  operate  for 
approximately  520  seconds  total  firing  duration.  The  require¬ 
ment  for  55  missions  totaling  7.5  hours  of  cumulative  operating 
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FIGURE  1  --  SIMPLIFIED  ENGINE  CONTROL  SYSTEM ,  STAGED  COMBUSTION  ENGINE 
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time  with  varying  thrust  levels  represented  the  first  signi¬ 
ficant  requirement  for  a  reusable  liquid  rocket  engine.  The 
SSME  is  a  very  efficient  engine  with  a  specific  impulse  of 
approximately  455  seconds  at  rated  power  level  (FPL)  (470,000 
pounds  of  thrust  altitude).  The  two  stage  combustion  is 
approximately  99.96$  efficient.  The  main  chamber  combus tic- 
pressure  is  approximately  3000  psia  at  rated  power  level.  ihe 
turbopumps  are  direct  drive,  i.e.,  they  have  no  gear  trains. 
Another  interesting  feature  of  the  SSME  is  that  it  ^usps  head 
pressure  to  supply  the  energy  to  start  the  engine.  T'/ere  are 
no  start  tanks,  turbine  spinners,  or  pyrotechnics  ir/olved’  in 
engine  start.  Electrical  igniters  are  supplied  in  tlyJ  fuel  and 
oxidizer  preburner  and  the  main  combustion  chamber^ to  initiate 
combustion.  / 

Each  of  the  rocket  engines  operates  at  a  fixture  ratio  of  6:1, 
and  over  a  throttling  range  of  109$-6jj£"nPL.  This  provides  a 
higher  thrust  level  during  lift  of-f'  and  the  initial  ascent 
phase,  and  allows  orbiter  acceleration  to  be  limited  to  3g 
during  the  final  ascent  phase.  The  engines  are  gimballed  to 
provide  pitch,  yaw,  and  roll  control  during  orbiter  boost  phase. 

SSME  HYBRID  SIMULATION  MODEL 

The  first  stage  in  designing  a  feedback  control  law  is  to 
derive  a  set  of  differential  equations  that  describe  the  system 
response  as  accurately  as  possible.  The  SSME  hybrid^  model 
developed  by  Rocketdyne  on  the  AD-10/PDP-11  computer  is  the 
simulation  model  used  in  this  study.  The  engine  system  analog 
simulation  model  is  constructed  on  a  component  basis.  Indivi¬ 
dual  turbines,  pumps,  heat  exchangers,  and  fluid  flow  passages 
are  characterized  by  equations  defining  component  variations, 
dynamic  relationships,  and  interface  requirements.  The  system 
simulation  includes  all  static  and  dynamic  formulations  that 
are  considered  of  importance  in  accurately  representing  overall 
start,  mainstage  control,  and  cutoff  behavior  of  the  engine. 

The  hybrid  model  contains  three  subsections  that  are  used  to 
develop  the  basic  building  blocks  to  represent  all  phenomena  of 
significance.  The  engine  fuel  flow  system  involves  all  physi¬ 
cal  processes  where  hydrogen  in  a  gas,  liquid,  two-phase,  or 
super-critical  state  is  handled.  The  oxidizer  flow  portion  of 
the  engine  system  involves  the  physical  processes  where  oxygen 
is  handled  by  the  engine  system  prior  to  being  involved-  in  any 
combustion  process.  The  hot-gas  portion  of  the  engine  system 
simulation  includes  all  component  processes  that  involve  hand¬ 
ling  oxygen/hydrogen  combustion  products.  This  model  contains 
simplifications  which  include  perfect  gas  law  assumptions 
instead  of  National  Bureau  of  Standards  Property  Tables,  curve 
fitting  of  some  functions,  and  linearization  of  second-order 
nonlinear  effects.  A  schematic  of  the  simulation  program  is 
shown  in  Figure  2. 
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The  analog  model  description  represents  SSME  dynamics  and  is 
applicable  for  simulating  engine  start,  throttling,  and  cutoff 
dynamics  at  nominal  conditions.  The  analog  model  differential 
equations  are  used  to  develop  the  control  linear  model  used  in 
subsequent  sections  of  the  paper. 

CONTROLLER  REQUIREMENTS 

The  SSME  control  system  is  designed  to  meet  or  exceed  all  per¬ 
formance  control  requirements  set  by  the  Contract  End  Item 
(CEI)  specification.  Typical  CEI  requirements  encompass  the 
performance  criteria  for  startup,  mainstage,  and  shutdown 
(Table  1).  As  can  be  observed  in  Table  1,  all  CEI  control 
requirements  are  met  or  exceeded  by  the  present  controller. 

For  this  study,  CEI  design  requirements  will  be  the  ultimate 

goal  with  the  added  goal  of  minimization  of  temperature  gradi¬ 
ents  (F/s)  in  the  preburners.  Excessive  temperature  gradients 

(spikes)  during  engine  start/cutoff  can  cause  localized  crack¬ 
ing  of  turbine  blades  (e.g.,  shortening  of  turbine  blade  life). 

In  summary,  the  SSME  controller  requirements  for  this  study  are: 

1.  Design  a  multivariable  control  loop  which  will 

be  capable  of  performing  control  functions  on  a 
simulated  (hybrid)  SSME  and  demonstrating  super¬ 
ior  performance  against  existing  SSME  controls. 

2.  Develop  the  above  multivariable  control  loop 
with  the  added  design  goal  of  minimizing  turbine 
blade  metal  thermal  gradients  (thermal  shock; 

F/s)  on  the  SSME. 

MODELING  TECHNIQUES 

Introduction 


The  solution  for  a  nonlinear  set  of  differential  equations  is 
difficult  or  impossible  to  solve  analytically.  An  approach 
that  will  allow  an  approximate  solution  to  the  process  equa¬ 
tions  taken  from  the  SSME  hybrid  model  is  called  small  distur¬ 
bance  -theory,  small  perturbation  theory,  or  lineari¬ 
zation.  Initially,  a  set  of  steady  state  operating  con¬ 

ditions  is  found  from  an  engine  balance.  These  determine  the 
values  of  the  state  and  control  variables  needed  to  maintain 
the  engine  at  the  steady  state  operating  conditions.  The  pro¬ 
cess  equations  are  then  linearized  about  the  engine  design 
point.  The  following  sections  describe  the  linearization  pro¬ 
cedure  for  a  general  set  of  nonlinear  equations  and  then  apply 
the  result  to  the  SSME  analog  equations  discussed  previously. 
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CEI  REQUIREMENTS  ROCKETDYNE  PERFORMANCE 
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Linearization  Procedure 


Since  there  are  eight  first-order  differential  equations  that 
are  necessary  to  describe  the  fuel  preburner  process  dynamics, 
the  linearized  system  will  be  eighth  order  with  two  control 
variables.  The  six  state  and  two  control  variables  are  shown 
in  Figure  3  with  definitions  given  in  Table  2. 

The  general  linearization  procedure  can  be  carried  out  in  basi¬ 
cally  two  steps  for  any  set  of  nonlinear  equations.  They  are: 

1.  The  set  of  equations  must  be  written  in  the  form 

f(x,  x,  u)  =  0 

by  moving  all  the  terms  to  the  left  of  the 
equality. 

2.  Expand  f  in  a  multivariable  Taylor  series  expan¬ 
sion  about  some  reference  trajectory-given  by 
f(0,x  ,u)  =  0.  This  expansion  is  valid  only  for 
small  changes  in  the  state  and  control  varia¬ 
bles,  x  and  u,  respectively.  All  second  order 
terms  are  considered  small  and  are  dropped 

Ax  =  A  Ax  +  B  Au 

where 


A  = 

B  = 


3f 

3x 


•||  (system  matrix) 


-  (control  ma 
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3x 
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SPACE  SHUTTLE  MAIN  ENGINE  SCHEMATIC  SHOWING  ENGINE  STATE,  CONTROL,  AND  AUXILIARY  VARIABLES 

FIGURE  3 


TABLE  2 


SUMMARY 

OF  SSME 
STATE 

PERFORMANCE  CONTROL  REQUIREMENTS 
,  CONTROL  AND  AUXILIARY 

VARIABLE  DEFINITIONS 

NO. 

STATE 

DEFINITION 

1 

DWFPO 

- 

Fuel  Preburner  Oxidizer  Flow  Rate 

(lbm/s) 

2 

SF2 

- 

High  Pressure  Fuel  Turbopump  Speed  (RPM) 

3 

PMFVD 

- 

Main  Fuel  Valve  Downstream  Pressure 
(PSI) 

4 

DWFPF 

- 

Fuel  Preburner  Fuel  Flow  Rate  (lbm/s) 

5 

DWFNBP 

- 

Primary  Fuel  Nozzle  Bypass  Flow  Rate 
(lbm/s) 

6 

PC 

- 

Main  Combustion  Pressure  (PSI) 

CONTROL 

7 

FPOV 

- 

Fuel  Preburner  Oxidizer  Valve  Setting 
(?) 

8 

CCV 

- 

Coolant  Control  Valve  Setting  (?) 

AUXILIARY 

9  TFP 


Fuel  Preburner  Temperature  (°R) 


and 
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where  |ref  means  evaluated  at  the  steady  state  operating  point. 
When  applied  to  the  sixth  order  model,  the  linear  model  becomes: 


A 


-2296.1  0 

53.5  -25.0 

0  0 

0  0 

21826.7  0 

0  0 


0 

0 

-24.4 

0 

0 

2. 


0 

-51.9 

0 

-983.4 

-21173.7 

0 


0 

0 

0 

0 

7684.9 

0 


0 

.081 

-32000. 

145.3 
0 

261.4 


B 


2 . 9xl05 
0 
0 
0 
0 
0 


0 

0 

0 

0 

0 

2.74x10 


4 


where 

AxT  =  ADWFP0,  ASF2,  APMFVD ,  ADWFPF,  APC ,  ADWFNBP 
AuT  =  AFP0V,  ACCV 

function.  For  example,  the  frequency  bandwidth  for  thrust  mod¬ 
ulation  will  be  significantly  different  than  that  for  turbine 
inlet  temperature  spikes. 
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Actuator  Dynamics 


For  most  purposes,  it  can  be  assumed  that  actuator  dynamics  do 
not  influence  the  system  or  plant  dynamics  significantly.  This 
is  equivalent  to  saying  that  the  actuator  dynamics  are  restric¬ 
ted  to  high  frequency  (large  eigenvalue)  regions  that  imply 
fast  response.  In  most  instances,  the  actuator  dynamics  are 
significantly  "faster”  than  the  plant  dynamics.  This  is  not 
true  for  the  SSME.  The  fuel  preburner  oxidizer  valve  and  the 
coolant  control  valve  have  open  loop  frequencies  of  -100  RAD/S, 
which  is  of  the  same  order  of  magnitude  as  other  system  dyna¬ 
mics.  Therefore,  the  actuator  dynamics  have  been  included  in 
the  overall  design  process.  The  eigenvalues  and  the  corres¬ 
ponding  eigenvectors  are  shown  in  Table  3«  As  can  be  seen,  the 
eigenvalue  corresponding  to  the  combustion  pressure  APc 
(7684.09  RAD/S)  is  positive,  which  implies  a  nonminimum  phase 
situation.  The  smallest  eigenvalue  (-24.9  RAD/S)  corresponds 
directly  to  ASF2,  the  fuel  turbopump  speed.  The  complex  conju¬ 
gate  pair  of  eigenvalues  corresponds  primarily  to  APMFVD,  the 
main  fuel  valve  discharge  pressure.  There  is  coupling 
between  APc,  A  DWFPF ,  ADWFNBP ,  and  APMFVD.  The  eigenvalue  that 
corresponds  to  DWFPO,  the  oxidizer  flow  rate,  is  coupled  to  the 
combustion  temperature.  The  two  actuator  eigenvalues  are  loca¬ 
ted  at  -100  RAD/S.  Mathematically,  the  control  problem  is  to 
move  the  nonminimum  phase  (  APc)  eigenvalue  to  the  left  of  the 
real  axis  without  triggering  undamped  oscillation  of  the  oxi¬ 
dizer  and  fuel  flow  rates,  which  in  turn  causes  the  fuel  pre¬ 
burner  temperature  to  oscillate.  The  actuator  dynamics  can  be 
approx iiuated  as  a  first-order  process  with  a  time  constant  of 
10  ms.  For  control  Au(t)  and  input  Ar(t),  the  differential 
equation  is 

Au ( t )  =  -  ~  u(t)  +  ^  Ar(t) 
where  t  is  the  time  constant. 

The  perturbation  model  for  the  fuel  preburner  oxidizer  valve 
(AFP0V)  and  the  coolant  control  valve  (ACCV)  is  as  follows: 


d 

"  AFPOV ( t ) " 

1  0  * 

’  AFPOV ( t )  ' 

dt 

ACCV  (t) 

0  -1 

ACCV(t) 

- 

°  T  J 

^  ■* 

AFP0VI  (t)* 
ACCV I ( t ) 
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where  AFPOVI(t)  and  ACCVI(t)  are  the  inputs  into^the  actuator. 
The  augmented  system  and  control  matrices  A  and  B  respectively 
are: 


A 

6x6 

0 

2x6 


B 

6x2 

C 

2  x  2 


B 


r  0  i 

6x2 

T 

2x2 


0  ' 
I. 

T 


C  =  -T 


and  0  is  the  null  matrix. 

The  linear  model  for  the  LQR  design  process  discussed  later  is 
therefore: 

Ax ( t )  =  A  Ax ( t )  +  BAu(t) 

where 

Ax1  =  [ADWFPO,  ASF2,  APMFVD ,  ADWFPF ,  APC, 

ADWFNBP,  AFPOV,  ACCV] 

Au*  =  [AFPOV I  ,  ACCV I ] 


Note  that  AFPOVI  and  ACCVI  are  perturbations  away  from  some 
nominal  operating  condition  for  the  actuator.  No  time  lag  is 
associated  with  the  fuel  preburner  temperature  and  is  consi¬ 
dered  to  be  a  nonlinear  function  of  the  oxidizer  fraction  in 
the  combustion  chamber  (ratio  of  oxidizer  flow  to  oxidizer  plus 
fuel  flow): 

TFP  =  T(FFP)  +  277. 4°R 

FFP  =  DWFPO  lbm/s 

DWFPO  +  DWFPF  lbm/s 

where 

DWFPO ( t )  =  ADWFPO ( t )  +  DWFPOref  lbm/s 
DWFPF(t)  =  ADWFPF (t)  +  DWFPFref  lbm/s 
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The  steady  state  operating  point  corresponds  to  the  following 
values  for  the  state  and  control  variables: 

DWFPOref  =  75.55  lbm/s 
SF2ref  =  3615.0  RPM 

PMFVDref  =  5963.0  psi 

DWFPFref  =  77.88  lbm/s 
Pcref  =  2995  psi 
DWFNBPref  =  65.48  lbm/s 
FPOVref  =  79.0* 

CCVref  =  99 . 9% 


OPTIMAL  CONTROL  TECHNIQUES 

INTRODUCTION 

Generally,  there  are  two  methods  of  designing  control  systems 
—  the  classical  and  the  modern  theory  approach.”  The  clas¬ 
sical  approach  usually  deals  with  single-input,  single-output 
linear  systems  in  the  frequency  and  Laplace  or  S-domains.  The 
modern  control  approach,  however,  deals  with  multi-input, 
multi-output  linear  systems  primarily  in  the  time  domain.  It 
is  the  latter  concept  that  is  used  herein. 

The  objective  of  modern  control  theory  is  concerned  with 
finding  a  suitable  control  law  usually  optimized  in  some  sense 
and  then  finding  a  controller  configuration  that  will  generate 
such  a  control  law.  The  control  law  is  not  constrained  to  take 
on  any  particular  form  but  for  most  purposes,  it  is  taken  so  as 
to  cause  control  deflections  proportional  to  some  error.  Typi¬ 
cally,  the  errors  of  interest  are  those  due  to  differences  in 
the  actual  value  of  variables  describing  the  process  such  as 
pressure  and  flow  rate  and  the  values  these  variables  take  for 
some  reference  or  equilibrium  conditions;  one  might  say  the 
controls  are  deflected  proportional  to  some  perturbation  away 
from  the  reference  condition.  The  problem,  therefore,  is  to 
select  gains  in  such  a  manner  that  multiplying  the  errors  by 
the  feedback  gains  provides  a  control  signal  in  such  a  manner 
as  to  maintain  desired  system  characteristics,  response  and 
stability.  The  method  used  to  determine  these  gains  is  based 
on  linear  system  theory  and  because  of  its  nature  is  called 
Linear  Quadratic  Gaussian  Theory  (LQG).^ 

The  LQG  design  procedure  is  basically  a  three  point  procedure: 

1.  Deterministic  Ideal  Response  Analysis  and  Design 

Step  1  involves  modeling  the  physical  situation  in  the  form  of 
a  set  of  mathematical  equations.  Usually  these  equations  are 
nonlinear  and  must  be  linearized  as  discussed  previously.  This 
model  assumes  no  uncertainty  in  the  modeling  of  the  plant 
(physical  situation)  or  measurements.  Exact  measurement  of  all 
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plant  state  and  measurement  variables  is  assumed.  Actuator  and 
plant  dynamics  are  also  assumed  to  be  known  exactly.  The  SSME 
deterministic  model  formulation  describes  the  behavior  of  the 
fuel  preburner  state  and  measurement  variables  as  well  as  the 
interaction  of  the  thrust  (combustion  pressure)  with  changes  in 
fuel  preburner  flow  rates.  The  preburner  temperature  taken 
here  as  an  auxiliary  variable,  is  assumed  to  be  a  nonlinear 
function  of  the  oxidizer  fraction  in  the  combustion  process  and 
the  fuel  temperature  in  the  lines.  For  design  purposes,  this 
nonlinear  equation  is  used  to  calculate  fuel  preburner  tempera¬ 
tures.  The  dynamics  of  the  temperature  are  of  such  high  fre¬ 
quency,  that  this  dynamic  effect  can  be  ignored.  The  measure¬ 
ment  or  output  variables  will  be  the  pressure  and  volumetric 
flow  rate  measured  at  the  fuel  flowmeter. 

2.  Stochastic  Estimation  Analysis  and  Design 

Step  2  introduces  uncertainty  into  the  linear  model  to  compen¬ 
sate  for  linearization  errors  as  well  as  errors  due  to  non¬ 
linear  model  uncertainty.  Most  plant  variables  in  a.  real 
system  cannot  be  directly  measured  by  conventional  techniques. 
Therefore,  a  state  filter  or  observery>,u  is  used.  Not  only 
are  plant  variables  assumed  uncertain  but  sensor  errors  are 
also  assumed  in  error.  These  errors  are  described  statis¬ 
tically  by  intensity  matrices  of  covariance  u . 

3.  Stochastic  Feedback  Control  System  Design 

From  Steps  1  and  2,  an  optimal  control  correction  from  the 
estimated  state  deviation  (error)  is  derived.  The  controller 
is  then  tested  in  the  linear  system  to  determine  how  good  the 
control  gains  are  in  meeting  design  goals. 

In  general,  a  control  system  is  designed,  for  each  set  of 
linearized  equations  that  describe  the  dynamics  of  the^SSME  in 
the  neighborhood  of  several  equilibrium  conditions.  For  this 
particular  problem,  the  mixture  ratio  of  6:1  at  rated  power 
level  (RPL)  is  of  interest  and  hence  is  used  as  the  reference 
condition  about  which  to  implement  the  controller  design..  The 
technique  for  designing  the  control  system  is  the  previously 
mentioned  LQG  theory,  which  provides  a  somewhat  systematic 
method  for  determining  a  set  of  feedback  gains.  In  most  cases, 
the  design  for  each  steady  state  condition  yields  a  set  of 
unique  feedback  gains.  The  resulting  designs  are  then  tested 
by  observing  the  time  responses  to  typical  input  commands.  In 
this  case,  these  commands  consist  of  a  step-type  combustion 
pressure  (thrust)  or  preburner  fuel  and  oxidizer  flow  rate 
changes . 
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THE  OPTIMAL  REGULATOR 


Introduction 


For  the  SSME  design  point  discussed  previously  dealing  with  the 
SSME  hybrid  model,  a  Linear  Quadratic  Regulator  (LQR)  design  is 
carried  out  that  provides  state  feedback  gains  to  the  control 
to  minimize  some  cost  function.'1  The  LQR  controller  always 
provides  a  control  law  that  will  drive  perturbations  or  errors 
from  the  ideal  state  and  control  variables  to  zero.  It  must  be 
understood  by  the  designer,  however,  that  the  LQR  formulation 
does  not  allow  constraints  to  be  put  on  the  control  or  state 
variables.  Such  constraints,  however,  can  be  applied  indirec¬ 
tly  by  the  designer  in  the  design  process.  In  general,  the  LQR 
design  procedure  selects  feedback  gains  in  such  a  manner  as  to 
minimize  the  "cost”  of  the  process  in  a  rigid  mathematical 
sense. 

In  the  following  sections,  the  general  formulation  of  the  LQR 
design  problem  is  presented.  In  addition,  a  procedure  will  be 
discussed  that  will  allow  the  designer  to  "constrain"  selected 
variables  in  order  to  stay  within  the  prescribed  limits  of  the 
engine. 

The  Regulator  Design  As  A  Tracking  Controller 

The  LQR  problem  centers  on  designing  a  controller  that  keeps 
errors  or  perturbations  small.  However,  the  real  problem  of 
interest  is  that  of  carrying  out  step  changes  in  the  main  com¬ 
bustion  pressure,  fuel  and  oxidizer  flow  rates  in  the  fuel  pre¬ 
burner  in  a  reasonably  short  time  while  retaining  good  dynamic 

characteristics  and  at  the  same  time  not  exceeding  specified 
engine  limits.  Although  these  two  problems  appear  different, 
they  can  be  cast  into  the  same  form  as  will  be  shown.  If  the 
engine  is  operating  at  one  steady  state  reference  point  and  it 
is  desired  to  change  to  another  point,  the  regulator  control 
can  be  used  by  simply  changing  the  steady  state  reference 
point.  Since  the  regulator  tends  to  drive  all  errors  to  zero, 
it  will  drive  the  engine  to  the  new  operating  point  with  new 
steady  state  conditions.  For  example,  suppose  it  is  desired  to 
change  the  thrust  level  in  15&  increments.  For  the  various 
thrust  levels,  there  will  be  a  unique  set  of  steady  state  oper¬ 
ating  conditions,  as  shown  in  Fi*gure  4.  This  new  set  of  steady 
state  operating  points  will  then  be  substituted  in  place  of  the 
old  reference  points.  The  change  in  the  variables  in  the  sys¬ 
tem  will  then  be  driven  to  zero,  and  hence  the  engine  will  be 
driven  to  the  new  operating  point  (new  thrust  level). 


41 


Reference  or  Desired  Variable  as  a  Function  of  Thrust 

Figure  4 


The  *  reference  map’  will  be  stored  in  the  computer  as  a  curve 
fit  or  table.  Note  that  thrust  is  assumed  as  the  independent 
variable  because  the  space  shuttle  guidance  computer  commands  a 
thrust  change.  It  is  possible  that  the  reference  variables 
could  be  a  function  of  other  independent  plant  variables  lik 
mixture  ratio,  for  example.  It  must  be  remembered,  however, 
that  the  ability  of  the  LQR  controller  to  accomplish  this  task 
is  a  function  of  how  well  the  linear  and  nonlinear  ®°del  cor¬ 
responds  at  the  steady  state  reference  point.  It  is  possioie 
that  a  n  change  in  the  reference  condition  may  exceed  the 
neighborhood  about  which  the  linear  perturbation  model  is 
valid.  In  that  case,  other  perturbation  models  oust  be 
designed  that  describe  the  motion.  This  'trim  map  methodo  gy 
of  designing  feedback  gains  for  various  operating  conditions  is 
well  known.  Of  course,  the  necessity  of  designing  new  linear 
models  depends  on  how  the  linear  and  nonlinear  transients  cor¬ 
respond,  given  the  same  initial  conditions.  It  is  intuitively 
obvious  that  'small'  changes  will  yield  good  correlation;  it  is 
also  true,  however,  that  large  changes  say  10*  of  the  steady 
state  values,  may  yield  results  that  do  not  represent  the  dyna¬ 
mics.  There  are  no  easy  methods  of  determining  the  linearized 
'neighborhood'  mentioned  above  other  than  nonlinear  simu¬ 
lation.  Hence,  although  a  regulator  design  is.  being  consi¬ 
dered,  a  tracking  controller  can  be  implemented  simply  by  vary- 
inc:  the  reference  control  in  some  prescribed  fashion  such  as  a 
step  change  used  here,  or  possibly  a  ramp,  sinusoid  or  some 
other  form. 
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Mathematical  Formulation 


1 1 


Briefly  stated,  we  wish  to  find  a  control  law  that  will  mini¬ 
mize  a  cost  or  penalty  function. 

j  =  1  r  ( AxTQ  Ax  +  A u T  R  Au)dt  (1) 
2  o 

where  Q>  0  and  R  >0  and  hence  drive  perturbations  in  x  and  u 
to  zero.  This  cost  function  is  subject  to  the  constraint  that 
the  variables  Ax  and  Au  must  satisfy  the  differential  equations 

Ax  *  A  Ax  +  B  Au  (2) 

where  A  and  B  are  constant  system  and  control  matrices, 
respectively.  A  quadratic  cost  functional  as  above  is  chosen 
because  it  penalizes  large  errors  much  more  severely  than  small 
errors.  As  is  clearly  seen,  the  first  term  in  Equation  (1) 
penalizes  errors  in  the  state  and  the  second  term  deviations 
away  from  the  reference  conditions  for  the  control.  There  are 
two  main  reasons  why  Equation  (1)  is  used.  They  are  that  it  is 
mathematically  tractable;  that  is,  the  theory  has  been  well 
established,  and  also  that  it  leads  to  optimal  linear  feedback 
systems. 

The  solution  to  the  problem  of  minimizing  Equation  (1)  subject 
to  Equation  (2)  is  given  by  a  constant  linear  feedback  matrix  K 
so  that 


Au ( t )  =  -K  Ax ( t )  (3) 

where  K  =  R“^B^S  and  S  >  0  satisfies  the  steady  state  matrix 
Riccati  equation. 

-SA  -ATS  +  SBR"1  BTS  -  Q  =  0 

Existence  of  a  unique  positive  definite  solution  to  the  LQR 
control  problem  is  guaranteed  provided  the  nxmn  "control¬ 
lability  matrix" 

[B:AB: - :  An"^B ] 

has  rank  n  where  n  is  the  dimension  of  the  A  matrix  and  m  is 
the  number  of  controls  (columns  in  the  B  matrix). 

Procedure 


In  general,  the  procedure  for  using  the  regulator  method  can  be 
divided  into  three  main  steps: 

1.  The  Q  and  R  weighting  matrices  must  be  selected 
for  the  cost  functional.  In  Equation  (1),  note 
that  Q  must  be  a  positive  semi-definite  matrix 
while  R  must  be  a  positive  definite  matrix. 
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The  gain,  matrix  K  must  be  determined  by  solving 
the  matrix  Riccati  equation  (4)  for  the  matrix  S 
and  substituting  the  result  in  Equation  (3). 

The  resulting  control  law  must  then  be  incor¬ 
porated  into  the  original  system  and  tested.  In 
general,  this  step  is  carried  out  in  two  sub¬ 
steps.  Initially,  the  controller  design  is  tried 
in  the  linearized  system  for  which  it  is 
designed,  and  if  satisfactory  characteristics  are 
achieved,  it  is  then  tried  in  the  original  non¬ 
linear  system.  Most  often,  this  controller  veri¬ 
fication  step  is  carried  out  using  digital  simu¬ 
lations  of  the  linear  and  nonlinear  systems  on  a 
high  speed  computer. 


Although  appearing  straightforward,  this  procedure  is  not  with¬ 
out  some  difficulties.  In  particular,  problems  can  be  encoun¬ 
tered  selecting  the  elements  of  the  weighting  matrices  in  the 
cost  function  Q  and  R.  Generally,  only  the  diagonal  elements 
are  used  with  first  guesses  equal  to  approximately  1/e^, 
is  .the  maximum  desired  error  in  the  variable  associated 
with,  that  term.  Here  then,  is  where  the  designer  can  initiate 
considerations  for  constraints  on  the  state  and  control  varia¬ 
bles;  that  is,  they  determine  a  "ball  park"  number  for  e.  Most 
likely,  however,  step  3  reveals  a  dynamic  response  of  the  con¬ 
trolled  or  closed-loop  system  that  is  not  that  desired,  and 
adjustments  of  the  weights  are  necessary.  Oftentimes,  in  order 
to  avoid  unnecessary  computer  simulations  of  Step  3  as  an  aid 

°/.Q  and  5  weiShting  matrices,  the  eigenvalues 
rrom  the  closed-loop  system  are  used  to  determine  various  per- 

formance  characteristics  such  as  rise  time,  time-to-half  ampli¬ 
tude,  etc.  If  these  characteristics  are  favorable,  that  is,  if 
the  closed-loop  system  response  is  "fast",  then  the  system 
response  is  simulated  on  a  computer.  If  the  closed-loop  system 
is  "slow",  then  the  weights  are  changed  to  achieve  a  faster 
response.  .  The  fast  and  slow  designations  are  related  to  the 
relative  time  it  takes  to  achieve  a  steady  state  condition. 


A  method  that  proves  to  be  successful  in  choosing  the  Q  and  R 
weighting  matrices  is  to  change  the  diagonal  terms  one  at  a 
time  and  note  changes  in  the  closed-loop  system  characteristics 
(eigenvalues).  Each  new  system  is  simulated  on  a  computer  to 
give  an  indication  of  how  changing  one  element  at  a  time  in  the 
Q  and  R  matrices  changes  the  closed-loop  time  responses. 

As  mentioned  earlier,  the  vehicle  limits  are  not  part  of  the 
LQR  problem.  They  are,  however,  part  of  the  overall  design 
process.  Because  this  system  will  be  used  as  a  tracking  con¬ 
troller,  some  quick  checks  can  be  made  to  determine  if  the  Con¬ 
trols  will  exceed  their  prescribed  limits  for  a  given  design. 
Since,  for  a  step  change  in  main  combustion  pressure  (flow 
rates),  all  the  states  except  for  the  biased  variable  will  be 
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zero,  the  initial  control  deflections  can  be  determined  by  mul¬ 
tiplying  the  feedback  gain  corresponding  to  the  pressure  change 
required  by  the  desired  change.  The  resulting  control  deflec¬ 
tion  should  be  within  the  vehicle  control  limits.  If  the  con¬ 
trol  limit  is  exceeded,  then  the  element  corresponding  to  the 
control  in  the  R  matrix  is  increased.  This,  in  effect,  pena¬ 
lizes  the  control. 

In  general,  increasing  diagonal  elements  in  the  Q  and  R  ma¬ 
trices  tends  to  penalize  the  variables  corresponding  to  those 
elements.  Therefore,  in  order  to  limit  a  state  variable, 
increase  the  diagonal  element  in  the  Q  matrix  corresponding  to 
that  state  variable.  In  order  to  limit  a  control  variable, 
increase  the  diagonal  element  in  the  R  matrix  corresponding  to 
that  controlled  variable. 

This  trial  and  error  approach,  augmented  by  the  knowledge  of 
which  errors  are  not  important,  all  come  into  play  when  selec¬ 
ting  weighting  matrices.  In  any  case,  a  regulator  solution  and 
hence,  the  corresponding  feedback  matrix  yielding  desirable 
closed-loop  system  dynamics  can  be  obtained  with  patience. 


CONTROLLER  DESIGN 


INTRODUCTION 


The  Linear  Quadratic  Gaussian  control  scheme  as  discussed  pre¬ 
viously  is  a  relatively  simple  controller  to  implement  in  a 
computer  control  system  once  the  linearized  models  are 
obtained.  For  this  paper,  only  the  LQ  regulator  design  metho¬ 
dology  and  implementation  will  be  discussed.  The  addition  of 
the  Kalman-Bucy  filter  is  not  difficult  and  the  overall 
"marriage’'  of  the 
filter  can  be  found 


optimal  state  feedback  with  the  Kalman-Bucy 
in  any  text  on  linear  optimal  control.  2 


The  design  point  of  10056  rated  power  level  with  the  reference 
conditions  for  the  six  states  and  two  controls  is  used  as  a 
design  point  to  determine  a  feedback  matrix,  not  necessarily 
the  only  one,  that  will  produce  the  desired  behavior  in  the 
neighborhood  of  the  design  point.  Ideally,  one  would  like  all 
the  feedback  matrices  found  to  be  the  same  so  that  a  constant 
gain  system  results.  In  most  cases,  however,  the  feedback 
matrices  are  different.  Consequently,  either  some  scheme  must 
be  incorporated  into  the  controller  that  senses  the  engine's 
state  and  condition  or  some  alternate  scheme  must  be  deve¬ 
loped.  Although  there  are  several  variations  of  gain  sche¬ 
duling  techniques,  most  require  some  method  of  estimating 
states  and  parameters  (Kalman-Bucy  filter). 

The  LQ  regulator  gains  are  found  from  the  commanded  thrust, 
which  is  commanded  by  the  space  shuttle  guidance  computer. 
From  simulation  studies,  a  'trim  map*  of  state  and  control 
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variables  as  a  function  of  thrust  level  can  be  used  to  ’trim* 
or  drive  the  system  to  the  new  thrust  level.  For  different 
thrust  levels,  there  will  be  a  unique  set  of  reference  values 
for  the  state  and  control  variables.  Theoretically,  if  the 
full  set  of  actual  state  and  control  values  are  at  a  reference 
condition,  the  thrust  level  should  be  unique.  Therefore,  if 
the  computer  commanded  a  1$  change  in  thrust  level  from  say  65$ 
RPL,  the  computer  would  search  a  set  of  steady  state  or  refer¬ 
ence  state  and  control  variables  that  are  curve  fitted  or  in 
table  form  for  the  commanded  or  desired  values  of  these  para¬ 
meters  needed  to  drive  the  system  to  the  new  commanded  thrust 
level.  The  feedback  gains  will  then  drive  the  system  to  the 
new  state  and  control,  and  hence  to  the  commanded  thrust  level. 

The  LQ  gains  are  found  from  a  curve-fit  or  table  lookup  of  the 
gains  versus  the  thrust  level.  For  small  changes  in  thrust 
level,  the  gains  should  vary  smoothly.  Although  there  is  no 
guarantee  that  this  will  occur,  it  has  been  found  for  vehicle 
control  that  LQ  gains  will  vary  in  a  smooth  manner  if  the  state 
and  control  weighting  matrices  do  not  need  large  adjustments 
for  good  control. ^ 3 

The  LQ  gains  determined  for  the  design  point  model  discussed  in 
this  paper  were  found  using  ORACLS.'1*  This  package  is 
well-known  among  the  modern  control  community  and  will  not  be 
discussed.  The  use  of  this  package  of  routines,  however, .  is 
recommended  because  of  its  simplicity  in  accessing  various 
matrix  manipulation  and  LQG  design  algorithms. 

The  following  sections  discuss  the  controller  design  for  the 
fuel  preburner  linear  system  derived  previously.  A  section 
discussing  controller  implementation  as  well  as  linear  simula¬ 
tion  results  is  included.  Also,  the  effect  of  limiting  the 
control  in  the  linear  simulations  is  discussed. 

LINEAR  MODEL  SELECTION 

From  past  studies,  it  was  determined  that  fourteen  states  and 
five  controls  present  in  the  SSME  could  be  used  to  develop  a 
linear  model  to  describe  the  engine  throughout  the  flight  enve¬ 
lope.  However,  since  the  control  objective  was  to  control  pre¬ 
burner  temperatures,  it  was  decided  to  concentrate  on  the  fuel 
preburner  since  thermal  spikes  there  are  more  severe.  Further 
modeling  experience  indicated  that  only  nine  state  and  three 
control  variables  could  have  any  dynamic  effect  on  fuel  pre¬ 
burner  temperature.  From  an  eigenstudy  of  the  open-loop  dyna¬ 
mics  of  the  system,  it  was  found  that  certain  states  had  vir¬ 
tually  no  effect  on  preburner  temperature  because  of  the  high 
frequency  character  of  their  associated  eigenvalues.  It  was 
also  discovered  that  the  associated  eigenvector  was  dominated 
by  that  state.  Intuitively,  this  means  that  the  particular 
state  at  a  high  frequency  equilibrates  long  before  other  varia¬ 
bles  in  the  system.  (This  type  of  analysis  can  be  thought  of 
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as  modal  decomposition  since  the  system  matrix  yielded  no 
repeated  eigenvalues  and  hence  the  Jordon  matrix*^  is  just 
the  matrix  of  eigenvalues.) 

After  eliminating  the  row  and  column  in  the  system  matrix,  the 
eigenvalues  and  eigenvectors  were  calculated  again.  In  the 
variables  deleted  from  the  linear  model  representation,  it  was 
found  that  all  high  frequency  dynamics  were  eliminated.  It  was 
also  discovered  that  the  states  needed  for  the  fuel  preburner 
temperatures  did  not  have  the  same  frequency  level  as  the  dele¬ 
ted  states.  The  final  linear  system  contained  six  state  and 
two  control  variables.  It  was  discovered  later,  however,  that 
the  actuator  dynamics  associated  with  the  two  controls  were  in 
the  same  frequency  range  as  the  plant  dynamics.  Since  there  is 
no  direct  coupling  between  the  actuator  equations  and  the  plant 
equations,  the  actuator  poles  were  simply  negative  of  the 
inverse  of  the  actuator  time  constant  (e.g.,  -  —  )  taken  as 
10  ms.  The  resulting  eighth  order  system  was  used  Tfor  the  con¬ 
troller  design  model. 

As  previously  discussed,  the  effectiveness  of  the  control  law 
depends  primarily  on  how  well  the  linear  system  represents  the 
dynamics  of  the  nonlinear  system,  since  it  is  almost  certain 
that  linear  systems  found  from  nonlinear  simulations  will  be 
inaccurate  due  mainly  to  the  fact  that  the  nonlinear  simula¬ 
tions  cannot  describe  all  dynamic  variations.  Also,  since  the 
linear  models  are  not  of  full-order,  i.e.,  all  the  state  and 
controls  possible  are  not  used,  the  linear  control  law  derived 
will  yield  some  steady-state  error.  Although  this  error  should 
not  be  severe,  it  might  be  advisable  to  include  elements  of 
further  compensation  in  the  controller.  The  easiest  method  to 
implement  is  to  implement  trim  integrators  on  those  states 
exhibiting  steady  state  error.  This  amounts  to  integrating  the 
steady  state  error  over  time  and  multiplying  that  by  an  inte¬ 
gral  gain. 

SINGLE  POINT  LQR  DESIGN 

Since  this  study  is  a  feasibility  study  into  the  applicability 
of  the  LQG  technique  toward  preburner  temperature  control,  only 
one  reference  or  design  point  was  chosen  about  which  to  derive 
a  linear  model.  The  typical  design  strategy  for  this  design 
point  linear  system  is  to  choose  a  set  of  Q  and  R  weighting 
matrices.  Realistically,  a  real  controller  must  incorporate 
gains  found  for  a  ’trim  map’  of  the  engine  operating  range 
(e.g.,  varying  thrust  levels).  For  design  points  that  are 
relatively  ’close’  to  one  another,  a  ’middle’  design  point  can 
be  used  to  determine  a  set  of  Q  and  R  weighting  matrices  that 
can  be  used  over  the  entire  set  of  operating  points  to  design 
control  laws.  The  LQR  design  for  the  one  design  point  was  car¬ 
ried  out  using  the  method  described  previously . The  resulting 
optimal  design  shows  that  a  step  change  in  Pc  (main  combustion 
pressure),  which  corresponds  to  a  change  in  thrust  level  with 
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all  other  deviations  zero,  resulted  in  good  dynamic  behavior 
especially  as  the  model  dynamics  influence  the  fuel  preburner 
temperature.  The  linear  system  simulation  run  for  a  Pc  chance 
from  nominal  of  -10  psi  for  a  time  interval  of  0  to  .1  s  is 
shown  in  Figure  5.  For  Pc  changes  greater  than  “10  Psi» 
tiple  step  changes  could  be  made  until  the  required  step  change 
is  made.  From  iterations  on  the  Q  and  R  diagonal  weighting 
matrices,  a  final  set  of  gains  was  found  using  the  weights 


Q  =  1.0,  1.0-107,  1.0,  1.0,  1.0-105  ,  1.0- 103  , 
R  =  1.0-106,  1.0- 103 


The  eigenvalues  from  the  closed-loop  system  show  considerable 
improvement  over  those  of  the  open-loop,  as  can  be  seen  below: 

EIGENVALUES 

OPEN-LOOP  (s-1)  CLOSED-LOOP  (s_1) 


The  gains  are: 
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The  improvement  is  mainly  in  speed  of  response  since  the  eigen¬ 
values  have  been  moved  to  higher  frequencies.  The  positive 
open-loop  pole  (7684.9  s"1)  has  been  eliminated  in  the 
closed-loop  system.  The  closed-loop  system  eigenvalues  show 
complex  conjugate  pairs.  These  eigenvalues,  however,  appear  to 
be  well-damped,  as  can  be  seen  in  the  simulation  results.  One 
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FUEL  PREBURNER  OXIDIZER  FLOW  RATE  VS  TIME 
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FIGURE  5  --  LINEAR  SIMULATION  RESULTS 
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FIGURE  5  --  LINEAR  SIMULATION  RESULTS  (CONTINUED) 
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FIGURE  5  --  LINEAR  SIMULATION  RESULTS  (CONTINUED) 
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FIGURE  5  --  LINEAR  SIMULATION  RESULtS  (CONTINUED) 
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FIGURE  5  --  LINEAR  SIMULATION  RESULTS  (CONTINUED) 
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FIGURE  5  --  LINEAR  SIMULATION  RESULTS  (CONTINUED) 


FUEL  PREBURNER  OXIDIZER  VRLVE  DYNAMICS  VS  TIME 
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FIGURE  5  --  LINEAR  SIMULATION  RESULTS  (CONTINUED) 


COOLANT  CONTROL  VALVE  DYNAMICS  VS  TIME 
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FIGURE  5  --  LINEAR  SIMULATION  RESULTS  (CONTINUED) 
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FIGURE  5  --  LINEAR  SIMULATION  RESULTS  (CONTINUED) 


COOLANT  CONTROL  VALVE  SETTING  7.  VS  TIME 
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FIGURE  5  --  LINEAR  SIMULATION  RESULTS  (CONTINUED) 
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FIGURE  5  --  LINEAR  SIMULATION  RESULTS  (CONTINUED) 
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FIGURE  5  --  LINEAR  SIMULATION  RESULTS  (CONTINUED) 


of  the  major  goals  in  this  study  is  to  control  fuel  preburner 
temperature  through  oxidizer  and  fuel  flow  rate  control.  Since 
fuel  flow  is  difficult  to  regulate  at  this  operating  point 
(i.e.,  the  coolant  control  valve  is  99.9%  open,  which  implies 
maximum  fuel  flow  rate),  fuel  control  will  yield  minimal  effect 
on  preburner  temperature. 

The  actuator  dynamics  significantly  affect  closed-loop 
response,  since  commands  from  the  feedback  control  law  do  not 
have  enough  time  to  reach  their  steady  state  values  before  the 
next  control  command  is  input.  A  relatively  large  commanded 
change  in  valve  position  will  not  be  seen  in  the  actuator  dyna¬ 
mics  response,  since  there  is  not  enough  time  for  the  valve  to 
respond  (e.g.,  10  ms  time  constant).  This  effect  can  be  seen 
directly  by  comparing  the  graphs  of  the  fuel  preburner  oxidizer 
and  coolant  control  valves  setting  commands  with  the 
corresponding  actuator  response.  For  large  perturbations,  a 
nonlinearity  due  to  saturation  of  the  valve  limits  the  amount 
of  valve  movement  that  can  be  accomplished. 

The  maximum  oxidizer  flow  rate  amplitude  is  about  10  lbm/s, 
which  shows  up  as  a  temperature  increase  of  about  200°R. 
However,  since  the  oxidizer  flow  rate  changes  damp  out  in  about 
.03s,  the  preburner  temperature  also  damps  out  as  fast. 

As  the  coolant  control  valve  position  is  increased  from  99.9%, 
the  main  fuel  valve  discharge  pressure  increases.  However,  the 
CCV  limit  of  100%  full  open  causes  saturation  of  the  valve  set¬ 
ting,  which  causes  oscillation  in  the  fuel  flow  rate.  But, 
since  the  amplitude  of  the  oscillation  is  about  0.5  lbm/s,  the 
effect  on  fuel  preburner  temperature  is  minimal.  A  small 
’ripple*  can  be  seen  in  the  temperature  transient  in  the  first 
10  ms.  The  commanded  change  in  main  combustion  temperature  of 
-10  psi  is  damped  out  in  about  5  ms,  although  an  overshoot  of 
2  psi  can  be  seen.  The  other  transients  associated  with  fuel 
flow  rate,  DWFNBP,  DWFPF  and  SF2  damp  out  in  2  ms,  the  mixture 
ratio  transient  has  settled  in  less  than  2  ms. 

Controller  Implementation  in  the  Nonlinear  System 

The  feedback  gains  found  earlier  are  designed  using  a  linear 
system  with  constant  coefficients.  The  variables  which  are 
involved  with  the  linear  system  are  error  or  perturbation  vari¬ 
ables.  Hence,  the  state  error  is  multiplied  by  the  gain  matrix 
and  the  perturbation  control  is  computed.  In  order  to  imple¬ 
ment  this  type  of  controller  in  an  actual  engine  or  in  this 
case  in  the  nonlinear  engine  model  simulation,  it  is  necessary 
to  convert  actual  variables  into  error  or  perturbed  variables, 
apply  the  controller,  then  convert  the  perturbed  control  com¬ 
mands  to  actual  control  commands. 

The  feedback  gains  are  implemented  as  in  Figure  6.  From  the 
nonlinear  system,  the  state  variables  x(t)  are  found.  The 
reference  or  steady  state  values,  xn(t),  are  then  subtracted 
from  the  current  values  to  give  'the  perturbed  values  Ax(t). 
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FIGURE 
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--  FLOW  CHART  SHOWING  LQR  CONTROLLER  IMPLEMENTATION 


The  control  law  is  then  used  to  generate  the  perturbed  control 
variables  Au(t).  This  value  is  then  added  to  the  reference 
control  values,  un(t),  to  give  the  actual  control  variable 
u(t)  which  is  then  limited  to  set  actuator  limits  of  100/5.  The 
feedback  matrix,  K,  is  the  LQR  feedback  matrix  described 
earlier. 

CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FURTHER  STUDY 

The  following  conclusions  can  be  drawn  regarding  the  feasi¬ 
bility  of  implementation  of  an  LQR  tracking  control  on  the  SSME. 

1.  The  feasibility  study  was  conducted  with  the 

assumption  of  perfect  state  knowledge.  Measure¬ 
ments  of  some  variables  may  not  be  available.  It 
is  possible,  however,  to  measure  enough  para¬ 

meters  so  that  the  system  will  be  completely 
observable  and  hence  a  Kalman  Filter  or  stochas¬ 
tic  observer  can  be  used  to  estimate  system 
parameters. 

2.  Linear  models  need  to  be  designed  at  different 

operating  points  before  any  conclusions  can  be 
made  as  to  the  control  of  preburner  temperature, 
mixture  ratio  and  thrust  for  different  engine 
operating  conditions. 

3.  A  linear  model  should  be  developed  which  includes 
both  the  oxidizer  and  fuel  preburners  as  well  as 
the  main  combustion  chamber  dynamics.  A  model  of 
this  sort  will  include  the  oxidizer  flow  dynamics 
which  from  this  study  affect  preburner  tempera¬ 
ture  significantly. 

This  study  has  shown  that  preburner  temperature  control  can  be 
accomplished  using  a  simplified  model  of  the  preburner  and  the 
Linear  Quadratic  Regulator  theory  assuming  full-state  availa¬ 
bility  for  a  typical  main  combustion  pressure  change.  Further 
research  in  this  control  problem  can  be  divided  into  the  fol¬ 
lowing  areas. 

1.  Low-order  linear  model  derivation  that  well 
represents  the  system  dynamics  at  various  opera¬ 
ting  points.  This  model  should  include  dynamics 
relevant  to  temperature  control  including  the 
fuel  and  oxidizer  preburners  and  turbines  as  well 
as  the  thrust  chamber  dynamics.  These  models  can 
be  derived  through  perturbation  studies  of  a  non¬ 
linear  simulation  model  or  by  some  piecewise  ele¬ 
ments  of  the  system  and  control  matrices  by  com¬ 
puter  differentiation  of  the  nonlinear  equation 
evaluated  at  various  operating  conditions.  Most 
models  generated  by  these  methods  are  too  large 


63 


and  do  not  contain  the  most  convenient  parameter¬ 
ization  of  the  dynamics.  Various  methods  can  be 
used  to  derive  this  linear  model  the  most  attrac¬ 
tive  of  which  involves  modal  decomposition  which 
simply  performs  an  eigenanalysis  of  the 
full-state  linearized  digital  SSME  dynamic 
model.  From  eigenvalue  and  eigenvector  studies, 
the  matrix  of  eigenvalues  should  yield  parti¬ 
tioned  submatrices  that  will  consist  of  high, 
middle,  and  low  frequency  dynamics.  The  high 
frequency  dynamics  are  assumed  to  be  equilibrated 
(Ax=0)  and  hence  reduce  the  model  order.  Some 
middle  frequencies  can  probably  be  eliminated 
from  further  study.  It  must  be  remembered,  how¬ 
ever,  that  the  partitioning  of  the  system  is 
dependent  on  the  control  designer’s  estimate  of 
the  frequency  range  of  the  control  function.  For 
example,  the  frequency  bandwidth  for  thrust  modu¬ 
lation  will  be  significantly  different  than  that 
for  turbine  inlet  temperature  dynamics.  Actuator 
dynamics  should  also  be  included  and  realistic 
time  delays  built  in  to  the  simulation  models  for 
control  activity. 

2.  The  set  of  low-order  linearized  models  are  then 
used  to  design  LQ-regulators . 

3.  A  sensitivity  study  of  the  system  should  then  be 
performed  to  yield  information  concerning  transi¬ 
ent  changes  due  to  varying  system  parameters  such 
as  thrust  for  example.  A  sensitivity  study 
should  also  yield  clues  to  the  type  of  gain  sche¬ 
duling  technique  to  be  implemented. 

4.  A  study  should  be  performed  that  will  develop  a 
set  of  Kalman-Bucy  filter  gains  (state  estimator) 
assuming  an  observable  system  for  each  linear 
model.  Gain  scheduling  techniques  should  also  be 
investigated. 

5.  The  Kalman-Bucy  filters  and  their  associated 
LQ-regulator  gains  should  then  be  tested  on  the 
linear  and  nonlinear  models  for  controller  effec¬ 
tiveness  over  the  flight  envelope. 

6.  A  feasibility  study  using  performance  results 
from  above  should  then  be  performed  to  determine 
any  hardware  incompatibilities  ignored  in  the 
above  steps.  If  the  control  designer  has  done 
his  job  correctly,  a  test  program  should  be 
investigated  to  test  the  controller  on  a  real 
engine . 
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This  sequence  is  not  meant  to  be  all  inclusive.  Any 
changes  in  this  sequence  will  be  a  function  of  experience 
as  the  control  study  progresses.  This  sequence  if  fol¬ 
lowed  should  yield  an  excellent  multivariable  control  law 
for  the  SSME  that  will  significantly  increase  engine 
performance. 
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ABSTRACT 

A  vehicle  suspension  dynamic  response  design  sensitivity 
analysis  and  optimization  technique  is  presented  and  illus¬ 
trated.  Dynamic  response  measures  included  in  the  formulation, 
for  use  as  the  objective  function  or  as  constraints,  include 
driver  absorbed  power,  driver  peak  acceleration,  and  suspension 
element  travel.  Design  parameters  available  to  the  designer 
include  suspension  spring  and  damper  characteristics,  suspen¬ 
sion  dimensions,  and  parameters  in  feedback  control  suspension 
subsystems.  An  adjoint  variable  technique  that  is  generally 
applicable  to  such  problems  is  employed  and  formulas  for  deriv¬ 
atives  of  dynamic  response  measures  with  respect  to  design 
parameters  are  derived.  Numerical  results  with  a  five  degree- 
of-freedom  vehicle  model  demonstrate  feasibility  of  the  method 
and  may  serve  as  a  guide  for  application  to  more  complex 
models. 

INTRODUCTION 

Vehicle  mobility  modelling  has  progressed  to  the  point  that 
performance  of  actual  or  concept  vehicles  can  be  analyzed  using 
proven  computer  simulations.  The  NATO  Reference  Mobility  Model 
(NRMM)  [1]  has  been  used  extensively  to  evaluate  vehicle  ride 
quality  and  overall  performance.  It  makes  use  of  driver 
absorbed  power  [2,3]  as  a  criteria  for  determining  maximum 
acceptable  sustained  vehicle  speed  over  terrain  and  driver 
acceleration  as  a  criteria  for  maximum  acceptable  vehicle  speed 
over  an  obstacle.  More  recently,  the  Dynamic  Analysis  and 
Design  System  (DADS)  methodology  has  been  used  in  detailed 
analysis  of  vehicle  dynamic  response  to  gun  firing,  passage 
over  terrain  and  obstacles,  and  weapon  station  stabilization 
[4,5]. 

To  date,  the  mobility  modelling  methods  noted  above  have  been 
used  almost  exclusively  for  evaluation  of  vehicles.  Little 
effort  has  been  devoted  to  extending  these  techniques  as  design 
tools.  Methods  of  control  system  optimization  [6]  have  been 
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developed  for  mechanical  system  optimization  [7],  but  have  not 
yet  been  brought  to  bear  on  vehicle  system  design  optimization. 
The  purpose  of  this  paper  is  to  present  and  illustrate  a  method 
that  can  be  used  to  build  upon  modelling  methodology  to  obtain 
a  design  tool. 

In  order  to  be  concrete  in  presentation  of  the  method,  a  five 
degree-of-freedom  vehicle  model  that  is  similar  to  idealiza¬ 
tions  used  in  the  NRMM  is  employed  as  an  illustrative  example. 
System  equations  are  derived  using  Lagranges  equations  to  yield 
second  order  differential  equations  that  depend  on  system, 
design  parameters.  These  equations  are  then. reduced  to  first 
order  form,  to  allow  easy  development  of  design  sensitivity 
equations.  Absorbed  power  calculations,  using  a  system  of 
first  order  differential  equations,  is  then  reviewed  as  a  key 
ingredient  in  design  problem  formulation. 

A  vehicle  design  optimization  formulation  is  presented  to 
minimize  driver  absorbed  power  on  a  nominal  road,  subject  to 
bounds  on  absorbed  power  on  a  rough  road,  driver  peak  accelera¬ 
tion  over  a  discrete  obstacle,  suspension  jounce  and  rebound 
travel,  wheel  hop,  and  limits  on  design  parameters.  An  adjoint 
variable  method  borrowed  from  optimal  control  theory  [6]  is 
then  derived  for  calculation  of  design  derivatives  of  vehicle 
dynamic  response  measures.  An  iterative  optimization  algorithm 
[8]  is  then  outlined  for  vehicle  design  optimization.  Its  use 
for  solution  of  the  five  degree-of-freedom  example  is  then 
illustrated. 


VEHICLE  MODEL  AND  EXCITATION 

The  vehicle  model  used  in  this  study  is  a  five  degree-of-free¬ 
dom  (5-DOF)  plane  linear  model  shown  in  Fig.  1 .  The  general¬ 
ized  vehicle  coordinates  are  the  passenger  seat  displacement 
x-j  ,  the  vehicle  body  vertical  displacement  X2  and  rotation  X3 , 
and  the  front  and  rear  wheel  vertical  displacements  X4  and  X5. 
The  spring  stiffnesses  are  denoted  by  k^  to  k5  and  damping 
coefficients  are  c-|  to  C5 .  Lengths  are  denoted  by  L-j  to  L4. 
The  functions  f-|  (t)  and  f2(t)  represent  displacements  of  the 
front  and  rear  wheels  due  to  undulation  in  the  road  surface  on 
which  the  vehicle  is  traveling. 

The  equations  of  motion  for  the  vehicle  are  derived  using 
Lagrange's  equations  of  motion.  In  matrix  form,  they  are 


[M]x  +  [C]x  +  [K]x  =  q(t) 


(1) 
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-k< 


"L4k1 


(k^ +k2+k3) 

(L^k1+L2k2-L3k3) 

-k2 

"k3 

[K]  = 

Symmetric 

(L^k1+L2k2+L3k3) 

"L2k2 

(k2+k4) 

L3k 

0 

q(t)  = 


0 

0 

0 


k^fjCt)  +  c^Ct) 
k^f2<t)  +  c^i^Ct) 
If  one  defines 


z.  =  x. 
x  1 


z5+i  “  xi 


i  =  1,2,3, 4, 5 


(k3+k5) 


(4) 


(5) 


(6) 


then  Eq.  1  and  initial  conditions  can  be  transformed  to  the 
first-order  form 


z(t)  =  f(t,  z,  b) 
z(0)  =  0 

where 


(7) 


z(t)  —  [ Zj  ,  z2 ,  •  •  •  z-\  n ] 


'10 


(8) 
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and 


f (t,z,b) 


b  —  [b^  »t“2>  •  •  —  [^-j  * ^2 >c-|  »c2,c3^  (9) 

z6 
z7 
z8 
z9 
Z1 0 

£j-(-b1  z1+b1  z2+L4b1  z3-b4z6+b4z7+L4b4zg) 

— (b1 z1 - (b1+b2+b3) z2- (L4b1+L2b2-L3b3) z3+b2z4+b3z5 

b4z6"^b4+b5+b6^z7"^L4b4+L2b5”L3b6^z8+b5z9+b6z10} 
l-{ L4b 1 z 1 - (L4b 1 +L2b  2- L3b  3 ) z  2- ( L4b 1 +L2b  2+L^b  3 ) z  3 

+L2b2z4~L3b3z3+L4b4Zg- (I^b^L^b^-I^b^Zy 

“(L4b4+L2b5+L3b6^z8+L2b5z9“L3b6z10} 
n^^b2z2^b2b2z3_b2z4~^4^ z4”^1 ) 
+b5z7+L2b6z8”b5z9”c4^  z9”^1 
mJ{b3z2”L3b3z3~b3z5"k5^z5_f2^t^  ) 
+b6z7‘L3b6z8"b6z10"c5^z10”f2(t^} 

(10) 
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Numerical  values  of  the  model  parameters  used  in  this  study  are 
given  in  Table  1 .  The  spring  and  damping  coefficients  given 
for  the  driver's  seat  and  suspension  will  be  the  starting 
values  for  optimization. 


Table 

1  . 

System 

Parameters 

Index 

1 

2 

3 

4 

5 

Generalized 
masses  m^ 

(lb  or  Id- in- sec  ) 

290 

4500 

41000 

96.6 

96.6 

Spring  coeff . 
ki 

(lb/in. ) 

100 

300 

300 

1500 

1500 

Damping  coeff. 
ci 

(lb-sec/in. ) 

10 

25 

25 

5 

5 

Distance,  L-; 

(in) 

120 

40 

80 

10 

— 

In  this  vehicle  model,  the  tires  are  modeled  as  point 
followers,  which  are  always  in  contact  with  the  ground.  If  a 
tensile  force  occurs  between  wheel  and  ground,  wheel  hop  would 
actually  occur.  In  the  design  formulation,  a  constraint  pre¬ 
cluding  wheel  hop  is  imposed. 

In  vehicle  dynamic  response  analysis,  the  input  road  conditions 
are  quite  important.  Dynamic  response  of  the  vehicle  depends 
strongly  on  the  vertical  displacement  history  of  the  wheels  on 
the  road  surface.  Typical  road  conditions  are  defined  as  a 
sinusoidal  undulation,  with  amplitude  yo  and  variable  half¬ 
wavelength  Aj_.  The  front  tire  displacement  v(y)  at  a  location 
y  is  defined  as 


i-1 


yQ[1  -  cos- - ^-]  ,  yi_1  <  y  <  y*",  i  is  odd 


v(y)  = 


yQ[1  +  cos^ 


i  (y-.y —1} ,  yi_1 

i 


<  y  <  y  ,  l  rs  even 


.  i 

where  y  is  a  coordinate  measured  along  the  road  and  y1  =  )>  A. . 

j=1  J 
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If  the  constant  speed  of  the  vehicle  is  denoted  by  s,  the 
elapsed  time  between  front  and  rear  tire  encounter  of  the  same 
point  on  the  road  surface  is  t0  =  L-|/s,  where  Li  is  the 
distance  between  front  and  rear  wheels.  Then, 

yg[1  -  cos  w^(t-t^  ^ )  ] ,  t^  ^  <  t  <  t'*", 
v(t)  =  . 

yQ[1  +  cos  w^t-t1'1  )  ]  ,  t1"  <  t  <  t1, 

where  w^  =  -ns/l±  and  t*-  =  yi/s.  The  vertical  displacement 
function  for  the  front  wheel  can  therefore  be  defined  as 


i  is  odd 

OD 

i  is  even 


v (t)  ,  0  <  t  <  T 

f n (t)  = 

0  ,  otherwise 


(12) 


where  Ti  is  the  time  at  which  the  road  undulation  ceases.  The 
vertical  displacement  of  the  rear  wheel  has  the  same  value  as 
that  of  the  front  wheel,  but  with  a  time  lag.  That  is, 

f2(t)  -  to  <  t  <  T,  +  to  0  3) 

In  this  study,  three  different  road  profiles  are  used  for 
determining  the  vertical  displacement  functions  of  the  tires, 
as  shown  in  Fig.  2.  Profile  (1)  is  a  continuous  sinusoidal 
curve  with  a  constant  half-wavelength  of  480  in.  and  an 
amplitude  of  yQ  =  2  in.  This  profile  represents  a  smooth  high¬ 
way  condition.  Profile  (2)  is  a  continuous  sinusoidal  curve 
with  a  constant  half-wavelength  of  60  in.  and  an  amplitude  of 
2  in. ,  which  represents  a  rough  highway  condition.  Profile  (3) 
is  a  combination  of  two  sinusoidal  curves  with  different  half¬ 
wavelengths  i.  i  =  360  in.  and  %2.  ~  144  in.  and  an  amplitude  of 
5  in.  This  profile  represents  a  severe  bump  condition. 

The  vehicle  speeds  for  each  road  profile  are  as  follows: 

960  in/sec  (54.5  mile/h)  for  profile  No.  1 

616  in/sec  (35.0  mile/h)  for  profile  No.  2 

450  in/sec  (25.6  mile/h)  for  profile  No.  3. 
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v(y) 

I 


(a)  profile  no.  1 


(c)  profile  no.  3 


Figure  2.  Road  surface  profiles 
ABSORBED  POWER 


Absorbed  power  is  a  measure  of  the  rate  at  which  vibrational 
energy  is  absorbed  by  a  human  and  is  a  quantity  that  may  be 
used  to  determine  human  tolerance  to  vibration  when  a  vehicle 
is  negotiating  rough  terrain.  The  "absorbed  power"  concept  was 
developed  in  the  1960's (2),  although  it  has  not  been  widely 
used  until  recent  years(l).  In  this  study,  the  absorbed  power 
in  the  time  domain  is  used  as  a  measure  of  ride  comfort. 


Absorbed  power  equations  given  in  the  NATO  Reference  Mobility 
Model  [1,2,3]  are  as  follows: 


where 


p(t)  =  g(p ,z,b) ,  0  <  t  <  x 

P(0)  -  0 

T 

p(t)  =  [ p i  ,P2  » • • • >Py] 


(14) 
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-29.8p1  -  497.49z1  -  100p2 


g  = 


1  Op-j 

1736. 9p1  -  108p4 

lOOp^  -  35.19P.J  -  39.1p4 

-315. 7p1  +  34.0956p4  +  171 .075p6 
-80.0p1  -  91.36p4  -  30.28p5 

Pi  -  0.108p4  +  0.25p^  -  6py 


(15) 


Here,  zi  is  the  vertical  acceleration  of  the  driver's  seat  of 
the  vehicle,  in  units  of  g's.  Thus,  from  Eq.  10, 


m]g{b1 (z1"z2"L4z3)  +  b4(z6'z7‘L4z8)} 


(16) 


where  g  is  the  gravitational  acceleration. 

The  average  absorbed  power  AP  is  defined,  in  units  of  Watts,  as 

AP  =  1  /  {p1(t)-0.108p4(t)+0.25p6(t)-p7(t)}2dt  (17) 

Using  the  solution  of  Eq.  14,  absorbed  power  can  be  calculated 
numerically. 


THE  OPTIMAL  DESIGN  PROBLEM 

With  equations  of  motion,  absorbed  power  equations,  and  terrain 
displacement  functions  defined  for  the  wheels,  one  can  now 
define  the  optimal  design  problem.  It  is  desired  to  make  the 
driver  as  comfortable  as  possible,  over  a  range  of  road  condi¬ 
tions  and  traveling  speeds.  Thus,  the  design  objective  is  to 
minimize  the  absorbed  power  of  the  driver  on  the  smooth  highway 
condition,  by  adjusting  suspension  parameters  of  the  vehicle, 
subject  bounds  on  the  following: 

i)  absorbed  power  and  wheel  hop  at  the  rear  wheel  on  the 
rough  highway  condition, 

ii)  maximum  vertical  acceleration  of  the  driver's  seat,  wheel 
hop  at  the  rear  wheel,  and  the  rattle  space  between  seat 
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and  body,  body  and  wheels,  and  wheels  and  ground  on  the 
bump  road  condition, 
iii)  design  variables. 

The  optimal  design  problem  may  be  written  as  follows: 


t  1 

Minimize  tn  =  /  G(p  (t))dt  (18) 

0  0 

where 

G(p)  *  {p1 (t)-0.108p4(t)+0.25p6(t)-p7(t)}2/i  (19) 

subject  to  the  state  equation  of  Eq.  7,  the  absorbed  power 
equation  of  Eq.  14,  and  the  constraints 


T 

’h  ■  / 
1  0 

G(p2(C))dt 

I 

-  61 

<  0 

(20) 

1|>2  =  2 

;5“£2^  "  0 2 

,  <  0, 

0  < 

t  <  T 

(21) 

♦3 = ' 

•(Z2-f2(t)}  - 

■  e3  < 

0,  0 

<  t 

< 

T 

(22) 

*4  = 

S7'b1<z2+V 

:3-Zi)+b4(z 

7+L4z 

,3. 

'8 

-  <  0,  0<t<T, 

(23) 

*5  = 

|  Z2"tL^Z2~  z  -j  | 

-  05 

<  0, 

0  < 

t 

<  T, 

(24) 

+6  = 

|  z^-z^-I^z^  | 

-  6  6 

<  0, 

0  < 

t 

<  x. 

(25) 

*7  = 

|  z^-Z2"L3z3  j 

'  07 

<  0, 

0  < 

t 

<  X, 

(26) 

*8  " 

l*4“fl(t)l  - 

6  8  < 

0  , 

0  < 

t 

<  X, 

(27) 

*9  =  1 

z\  -  f2(t)  - 

6  9  < 

0,  0 

<  t  1 

<  • 

r 

(28) 

*10  = 

“  (  z^~  f 2 ( t ) ) 

-  e10  <  0. 

0  < 

t 

<  T 

(29) 

where  6^,  i=  1,2,. ..,10  are  constraint  bounds  and  z3  and  pj 
are  state  and  power  solutions  and  j  is  the  road  profile  number, 
defined  as  j~1  for  the  smooth  highway,  j=2  for  the  rough  high¬ 
way,  and  j=3  for  the  bump,  and  the  constraints  of  design 
bounds 
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i-1,2,3,4,5,6 


(30) 


bb  <  b.  < 
1  i 


lU 

bi’ 


where 

9 i  =  maximum  allowable  absorbed  power. 

02,  9 9  =  maximum  allowable  upward  distance  between 
rear  wheel  and  road  surface.  Since  the 
tire  should  be  always  in  contact  with  the 
road,  this  value  may  be  static  deflection 
of  rear  tire,  neglecting  dynamic  effect  of 
tire. 

93»  910  “  maximum  allowable  downward  distance  be¬ 
tween  rear  wheel  and  road  surface. 

64  =  maximum  allowable  acceleration  at  driver's 

seat. 

65  =  maximum  allowable  distance  between  driver's 

seat  and  chassis. 

0g,  07  =  maximum  allowable  distance  between  body, 
and  front  wheel  and  rear  wheel,  respec¬ 
tively. 

0g  -  maximum  allowable  distance  between  front  wheel 
and  road. 

The  values  of  9j_  and  design  bounds  for  this  study  are  given  in 
Tables  2  and  3. 


Table  2.  0  ±  Values  Table  3.  Design  Bounds 


L 

U 

91 

3.5 

watts 

i 

bi 

bi 

0  9 

1.1127  in 

1 

251b/in 

5001b/in 

0  3 

3.0 

in 

1 OOOlb/in 

j 

2 

1 OOlb/in 

e4 

350 

in/sec 

3 

1 OOlb/in 

1 OOOlb/in 

05 

2.0 

in 

4 

1 .01b-sec/in 

501b-sec/in 

e6 

5.5 

in 

5 

2.51b-sec/in 

801b-sec/in 

9  7 

5.5 

in 

6 

2.51b-sec/in 

801b-sec/in 

e8 

2.0 

in 

09  1.1127  in 

010  3 . 0  in 
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DESIGN  SENSITIVITY  ANALYSIS 


Dynamic  systems  treated  here  are  described  by  a  design  variable 
vector  b  *=  [b-|  ,b? , . .  .bg]^  and  two  state  variable  vectors, 
z(t)  =  Ul  (t),z2(t)  , .  ..zio(t)]T  and  p(t)  =  ,  ,  J  , 

[pi (t) ,P2(t) , . . .p7(t) ]T,  which  are  the  solutions  of  initial 
value  problem  of  Eqs.  7  and  14,  rewritten  here  as 


z(t)  =  f(t,z,b),  0  <  t  <  t 

z(0)  =  0 

P(t)  =  g(p,z,b),  0  <  t  <  t 

p(0)  -  0 

The  design  problem  may  be  written  in  terms 
the  form 

i|/.  =  /  G(p1(t))dt  -  0.,  i=0,1 

1  0  1 

where  Bq  =  0  and 

^i  =  hi(ti,z,b)  -  0i,  i=2,3,. 


(31) 

(32) 

of  functionals  of 

(33) 

..10  (34) 


where  ti  is  determined  by  conditions 


!i(ti,z,b)  (jt  h^(t , z ,b)  J  ^ ^  0 , 


i=2,3, 


.10 

(35) 


The  dependence  on  design  variable  b  in  Eq.  33  arises  through 
the  absorbed  power  state  variable  p  =  p(t;  b) .  In  Eq.  34,  it 
arises  both  explicitly  and  through  the  displacement  state 
variable  z  =  z(t;  b) .  In  order  to  obtain  the  derivatives  of  ^ 
with  respect  to  b,  an  adjoint  variable  technique  [6,7]  may  be 
introduced.  For  the  i=0  and  1  , 


(36) 


where  subscript  denotes  differentiation.  One  may  introduce  an 
adjoint  variable  y  to  obtain  the  identity 
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(37) 


/  YT[p-g]dt  =  0 
0 

Taking  the  derivative  of  Eq.  37  with  respect  to  b,  one  has 

T  ,p 

-/  V  [Pb-8pPb-gzzb-gb]dt  =  0 
Integrating  the  first  term  by  parts  gives 

l  [i’+r'aplPb  +  vVb  +  vTgb  it 

-  YT(OpbCO  =  o  (38) 

since  pb(0)  =  0  from  Eq.  32. 

Since  Eq.  38  holds  for  any  function  y,  a  function  y  can  be 
chosen  to  satisfy  the  following  differential  equation  and 
terminal  condition: 

Y1  +  gpY1  =  Gp,  0  <  t  <  t 

,  i-1,2  (39) 

Yi(x)  =  0 

Equation  36  may  thus  be  written,  using  Eq.  38,  as 
di|> .  t  .  T  .  T 

3^  -  /  (-tt1  gfa-r1  g2^b)dt.  1-0.1  (40) 

Another  adjoint  variable  X  is  introduced  through  the  identity 

/  XT[z-f ]dt  =  0  (41) 

0 

Taking  the  derivative  of  Eq.  41  and  integrating  by  parts  gives 
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(42) 


/ 


0 


[xT+XTfz]zb  +  XTfb  dt  -  XT(t)zb(T) 


0 


since  z^(0)  =  0.  By  the  same  reasoning  as  above,  the 
additional  adjoint  equation  is  obtained  as 


X 1  +  f^X1  -  -gV,  0  <  t  <  T 


,  i=1 ,2 


(43) 


x  (t  )  -  0 


Thus,  Eq.  40  becomes 


d\|» .  t  .  T  .  T 

=  gb+X  fb^dt’ 


i=0 , 1 


(44) 


For  the  second  kind  of  functional,  Eq.  34  gives 


dip - 


3b1  *  hizIzb(ti)+z(ti)(ti)bl+hib+hittb  1=2,3  ”  ^5) 
Differentiating  Eq.  35  with  respect  to  b,  one  obtains 


[a.  -Hii  z(ti)](ti)h  +  A.  zb(tt)  +  a. 


(46) 


_ti 


Since  Eq.  35  is  to  determine  t^,  the  coefficient  of  (tj_)b 
cannot  be  zero  and  one  has 


^Pb  =  ‘  IT.  +n  .Z£(t  .)“  zb^ti^ 
ci  z 


n , 


fci 


_ b 

HZ  +HZ  f(ti) 

(47) 


Substituting  Eq.  47  into  Eq.  45,  one  has 


dip  - 


tv--  =  (h.  +3  .8.  }z,  (t.)  +  6  +  h.  ,  i*2 ,3  ,...10 

db  i_  r  d  1  1  Ik 


(48) 
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where 


f(t.)+h 

t 

+Q±  £(ti) 
.  Z 


i*2 , 3 , . . . , 1 0  (49) 


% 


Using  Eq.  42,  one  has  the  following  adjoint  problem 


x1  +  ty  -  o, 

X1(t1)  =  -hT  - 

z 


0  <  t  <  ti 

eini  H  Ht(ti#z,b) 
z 


i-2.3. 


.10 

(50) 


Then,  from  Eq.  44,  one  has 


d'I'j 

3^  -  hib  +  6l°lb 


(51) 


Computation  of  design  sensitivities  for  a  given  design  is  thus 
summarized  as  follows: 

For  i=0,1  (absorbed  power  functional): 

1.  Integrate  state  equations  of  Eqs.  7  and  14  from  t=0  to 
t,  where  t  is  the  time  when  dynamic  response  of  the 
vehicle  is  in  the  steady  state. 

2.  Integrate  adjoint  equations  Eqs.  39  and  43  backward  from 
t  to  0,  using  known  state  variables. 

3.  Calculate  design  derivatives  of  Eq.  44. 

For  i|i£,  i=2, 3 ,4 , . .  .1  0: 

1 .  Integrate  state  equation  of  Eq.  7  from  t=0  to  some  large 
time  for  given  road  condition  to  include  maximum  h^. 

2.  Find  t^  satisfying  Eq.  35  by  using  a  root  finding 
algorithm. 

3.  Integrate  the  adjoint  equation  of  Eq.  50  backward  from 
t=t£  to  0. 

4.  Calculate  design  derivatives  of  Eq.  51. 


The  cost  and  constraint  values  and  their  design  derivatives  for 
an  initial  design  b  =  [100,  200,  230,  2.8,  70,  22]T  are  cal¬ 
culated.  In  that  design,  only  ip-|  and  \p 5  are  violated. 

The  state  and  absorbed  power,  and  the  corresponding  adjoint 
equations  are  integrated  by  using  a  Runge-Kutta-Fiehlberg 
method  of  order  4  and  5,  with  a  relative  error  of  10“  and 
absolute  error  of  1 0“  .  The  error  tolerance  of  to  find  t2 

is  given  10“  . 
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To  check  the  validity  of  design  derivatives,  design  variations 
are  created  by  making  a  change  of  57«  in  each  design  variable. 
Results  for  i |>q,  iM  ,  and  ^5  are  given  in  Tables  4,  5,  and  6. 

The  predicted  variations  calculated  by  using  design  derivatives 
derived  here  show  very  good  correspondence  with  the  actual 
variation,  within  9%  difference. 

With  design  variation  fib  *  0. 05 [b-j  ,b2 , . .  .bg]^  and  multiples 
bi  =  b°  +  ifib,  i»1,2,3,4,  i.e.  ,  from  57«  to  207»  design  varia¬ 
tion,  the  actual  and  predicted  values  of  cost  and  constraints 
are  given  in  Figs.  3  to  13,  where  the  solid  line  is  actual 
value  and  the  dotted  line  is  the  predicted  value.  The  pre¬ 
dicted  values  show  good  correspondence  with  actual  values. 

ITERATIVE  OPTIMIZATION  ALGORITHM 

Many  mathematical  programming  algorithms  for  solving  optimiza¬ 
tion  problems  have  been  developed  and  evaluated  for  engineering 
design  optimization.  In  selecting  an  algorithm  to  be  used,  one 
may  base  his  choice  on  accuracy  of  the  result,  rate  of 
convergence,  time  required  for  computing,  computer  storage 
required,  and  compatibility  with  engineering  analysis  methods. 
It  is  not  possible,  however,  to  order  all  algorithms  in  one  and 
only  one  way  to  say  which  is  best.  An  algorithm  may  be  poor  as 
applied  to  a  certain  class  of  problems,  but  effective  in 
another  class.  Hence  it  may  be  necessary  to  have  a  reserve  of 
algorithms  and  to  choose  one,  depending  on  the  class  of 
problems  at  hand.  In  this  paper,  the  linearization  method  of 
Pshenichny  [8]  for  optimal  design  of  mechanical  systems  is 
used.  This  algorithm  was  originally  presented  in  the  Russian 
literature,  but  has  apparently  only  recently  come  to  the  atten¬ 
tion  of  workers  in  the  west.  Pshenichny  has  proved  convergence 
of  the  algorithm,  using  an  active-set  strategy  that  is 
essential  in  large  scale  mechanical  optimization  problems. 

The  general  mathematical  programming  problem  is  to  find  b  eRn 
to  minimize  fg(b),  with  constraints 


f^(b)  <  0,  i=1,2,...,m' 

f^(b)  =  0,  i=m'+1,...,m 


(52) 


where  f£(b),  i=0,1,...,m,  are  continuously  differentiable 
functions.  Note  that  f^(b)  =  0  is  equivalent  to  the  in¬ 
equalities  fi(b)  <  0,  and  -fi(b)  <  0.  Hence,  one  can  limit 
considerations  to  the  case  with  inequality  constraints.  Thus, 
one  wishes  to  minimize  fo(b),  subject  to  the  constraints 

fi(b)  <  0,  i=1  ,2 , . . .  ,m  (53) 
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Table  4.  Variation  of  cost  (i|>o)  with  57<>  change 
of  each  design  variable 


Design 

Variable 

Changed 

Current 

Cost 

Actual 

Change 

Predicted 

Change 

Percentage* 

Sensitivity 

b 

1 

6.568191 E—  0 1 

102.79 

b2 

6. 920365E-01 

1.507890E-02 

1  .4871 86E-02 

98.63 

b 

3 

6.935858E-01 

1 .66281 5E-02 

1 .686925E-02 

101 .45 

b 

4 

6. 708257E-01 

-6.1 31 887E-03 

-6.1 88627E-03 

100.93 

b5 

6.77961 OE-01 

1.003385E-03 

1.089884E-03 

108.62 

b 

6 

6. 745058E-01 

-2.451 777E-03 

-2.576824E-03 

105.10 

*percentage  sensitivity  is  defined  as  ^ "  x  100 


Table  5.  Variation  of  absorbed  power  constraint  (tl) 
with  57o  change  of  each  design  variable 


Design 

Variable 

Changed 

Current 

Constraint 

Value 

Actual 

Change 

Predicted 

Change 

Percentage 

Sensitivity 

b 

1 

2.79521 7E-01 

97.82 

b 

2 

6. 628705E-01 

3.448725E-03 

3.209578E-03 

93.07 

b 

3 

6.55861 5E-01 

-3.560305E-03 

-3.849868E-03 

108.13 

b 

4 

7.461 878E-01 

8. 676600E-02 

8.307421 E-02 

95.75 

b 

5 

8. 069681 E-01 

1 .475463E-01 

1 .490884E-01 

101.05 

b 

£ 

7.095821 E-01 

5.01 6029E-02 

4.964487E-02 

98.97 

Table  6.  Variation  of  rattle  space  constraint  between 
driver's  seat  and  chassis  (iK)  with  57o 
change  of  each  design  variable 


Design 

Variable 

Changed 

Current 

Constraint 

Actual 

Change 

Predicted 

Change 

Percentage 

Sensitivity 

5 

1 

b2 

2. 304623E-01 

1  .052389E-02 

1  .075973E-02 

102.24 

b3 

2. 265486E-01 

6.6101 85E-03 

6.89501 0E-03 

104.31 

b4 

2.05901 5E-01 

-1 .403692E-02 

-1 .41 6363E-02 

100.90 

b 

2.309891E-01 

1 .105070E-02 

1 .098409E-02 

99.40 

J 

b6 

2.1 50365E-01 

-4.901 91 6E-03 

-4.976799E-03 

101 .53 

Actual  Change 


o 


1 - 1 - 1 - 1 - 1 

0.0  S.O  10.0  15.0  20.0 

PERTURBATION (X) 


Figure  3.  Accuracy  test  of  design  sensitivity 
of  cost  ( \po) 
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Ar.tual  Chanae 


0.0  .  5.0  10.0  15.0  20.0 


PERTURBATION (X) 

Figure  4.  Accuracy  test  of  design  sensitivity  of  absorbed 
power  constraint  on  rough  highway  (<|cj) 


0.0  5.0  10.0  15.0  20.0 

PERTURBATION  [V.) 


Figure  5.  Accuracy  test  of  design  sensitivity  of  rear  wheel 
hop  constraint  on  rough  highway  (t2) 
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Figure  6.  Accuracy  test  of  design  sensitivity  of  rear  wheel 
J  penetration  constraint  on  rough  highway  ( 4* 3) 


:  Actual  Change 


Predicted  Change 


PERTURBATION  iV.) 

Figure  7.  Accuracy  test  of  design  sensitivity  of  peak 
acceleration  constraint  on  bump  road  (414) 
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:  Actual  Change 
:  Predicted  Change 


Figure  10.  Accuracy  test  of  design  sensitivity  of  rattle 
space  constraint  between  chassis  and  rear 
wheel  assembly  (^7) 


oj  - :  Actual  Change 


Figure  11.  Accuracy  test  of  design  sensitivity  of 
rattle  space  constraint  between  front 
wheel  assembly  and  road  (i|>g) 
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:  Actual  Change 


Figure  12.  Accuracy  test  of  design  sensitivity  of  rear 
wheel  hop  constraint  on  bump  road  (^9) 


:  Actual  Change 


Figure  13.  Accuracy  test  of  design  sensitivity  of  rear  wheel 
penetration  constraint  on  bump  road  (^iq) 
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BASIC  ASSUMPTIONS:  Let 

F(b)  -max  {0, (b) , . . . ,fm(b) }  (54) 

Note  that  F(b)  >  0  for  all  b  «  Rn.  Given  e  >  0,  define  the 
e-active  constraint  set 

A(b,e)  -  (isf^b)  >  F(b)  -  e,  i-1  ,2, . . .  ,m>  (55) 

(a)  Suppose  there  is  an  integer  N  y  0  such  that  the  set 

«  { b : 4>n (b)  <  *N(b®)}  (56) 

is  bounded,  where  b°  is  an  initial  design  and 
*N(b)  -  fo(b)  +  NF(b) . 

(b)  Suppose  gradients  of  the  functions  fj_(b),  i=0,1  ,2, . . .  ,m, 

satisfy  Lipschitz  conditions  in  there  exists 

L  >  0  such  that 

||f[(b’)  -  f^<b2) 1 |  <  Lllb1  -  b2||  (57) 

where  f[  -  [3f1/3b, , . . . ,3fi/9bn]T.  This  condition  is 

satisfied  if  fi  has  piecewise  continuous  first  deriva¬ 
tives.  .  ...  on 

(c)  Suppose  the  problem  of  quadratic  programming;  find  p  « 

to  minimize 

(fj(b).p)  +  \  IIpII2  <58> 

subject  to  the  linearized  constraints 

(f£(b),p)  +  ft(b)  <0,  i  «  A(b,e)  (59) 

is  solvable  with  any  b  *  and  there  are  Lagrange 
multipliers  u^(b),  i  *  A(b,e),  such  that 

V  u. (b)  <  N  (60) 

i«A(b,e)  1 
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THEORETICAL  ALGORITHM:  Under  the  above  hypotheses,  one  may 
state  the  following  theoretical  Linearization  Algorithm  of 
Pshenichny  [8]: 

Let  bO  be  an  initial  approximation  and  0  <  6  <  1 .  For  the 
kth  iteration, 

(TJ  Solve  the  quadratic  programming  problem  of  Eqs*  58  and 
59,  with  b  =  bk  and  solution  p*  =  p(b^) . 

(2)  Find  the  smallest  integer  i  such  that 


«(bk  +lT  pk)  <  *(bk)  -  6||pk||2 

21  21 

If  this  inequality  is  satisfied  with  i  »  ig,  let  =  2 


(61) 


-i 


0 


,  k+1  Kk  ,  k 

b  »  b  +  a,  p  . 


Under  the  basic  assumptions,  Pshenichny  [8]  has  proved  con¬ 
vergence  criteria  for  the  algorithm. 

NUMERICAL  ALGORITHM:  The  following  numerical  algorithm  is 
intended  for  solving  the  problem  of  minimizing  fo(b),  subject 
to  the  constraints  of  Eq.  52: 

Define 

F(b)  max{0,f1(b),...,ftn*(b),  |  fm'+1  (b)  |  , . . . ,  |  fm(b)  | } 

A(b, e)  =  (i:fi(b)  >  F(b)  -  e,  i  =  1,2,...,m’} 

B(b,e)  =  {i:  |fi(b) |  >  F(b)  -  e,  i  *  m'+1 , . • .  ,m} 

*N(b)  -  fo(b)  +  NF(b) 

Select  the  initial  approximation  b°,  Nq  sufficiently  large, 
eo  >  0,  and  0  <  6  <  1 . 

Step  1 .  In  the  kth  iteration,  solve  the  problem  of  finding  u 
to  minimize 


t(u)  -  J  |  |£g(bk)  +  I  Uiq(bk)||2-  I  U^iCbS 

i«A(bk,e)uB(bk,e)  i"A(bk,e)uB(bk,e) 

subject  to  u^  >  0,  i  e  A(b  ,e),  and  u^  arbitrary  for 

,  3 f .  3f.  T 

i  e  B(b  ,-e),  where  f|  =  [yg— , . . .  >y£~]  • 

1  n 


If  the  solution  u^  is  such  that  <j>  (u^)  =  then  set 

bk+1  =  b^,  efc+i  =  (1/2) e^,  and  Nfc+i  =Nk  and  return  to  Step  1. 

Otherwise,  let 
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(62) 


pk  -  -fj(bk)  -  l  ukq<bk) 

l«A(bk,e)uB(bk,e) 


and  go  to  Step  2. 
Step  2.  Set 


bk+1 


Kk  ,  k 
b  +  c*kp 


ek+1  "  ek 


where  <*k  is  chosen  equal  to 
for  which 


2q0 


and  qQ  is  the  smallest  integer 


V  (bk  +  3Pk)  ‘  4N,(bk>  '73  6  I  lpkl  I2 


2q 


2q 


Step  3.  If 


k  ,  v 

u..  +  l 


2  /  I 

ieA(b^,e)  ieBCb^.e) 


k  _  \  j  Ti  / 1_  k 


l“il\  >  >  l 


k  ,  v  I  k  | 

k  "  <■  .  Ui  +  i  .  lUll 

i«A(b  , e )  i«B(bK,e) 


Then  let  N^+l  =  Nfc.  Otherwise,  let 


j_i  -  l.  i  /  u.  H-  y  |u.  l 

k+1  I  L  v.  i  L  ,  1  L1 

i«A(bK , e )  i«B(b  , e ) 


-  2  /  l 


Step  4.  If  | |pk| |  is  sufficiently  small,  terminate.  Other¬ 
wise,  return  to  Step  1 . 

In  implementing  the  algorithm,  the  derivatives  f'  are  cal¬ 
culated  using  the  adjoint  variable  method,  which  requires  that 
the  state  and  adjoint  equations  be  solved  forward  and  backward 
in  time,  respectively.  This  is  a  substantial  amount  of  com¬ 
putation  that  must  be  carried  out  during  each  optimization 
iteration. 


* 


NUMERICAL  RESULTS 

Optimization  of  the  vehicle  suspension  is  carried  out  using  the 
linearization  method.  The  initial  design  shown  in  Table  7  is 
chosen  from  Ref.  7.  With  the  given  initial  design,  the  peak 
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acceleration  is  331.8  in/sec2  and  the  cost  is  0.847  watts. 

Also,  in  this  design,  constraints  on  absorbed  power  on  the 
rough  highway,  rear  wheel-hop,  and  rattle  space  constraint 
between  chassis  and  wheel  assemblies  are  violated.  This 
initial  design,  from  the  viewpoint  of  "ride  comfort"  and 
"safety,"  is  very  poor. 

After  the  7th  iteration,  the  rear  wheels  are  constantly  in 
touch  with  the  road  surface;  i.e. ,  the  wheel  hop  constraint  is 
satisfied.  Constraints  on  absorbed  power  and  rattle  space 
between  chassis  and  front  wheel  assembly  are  satisfied  after 
the  21st  iteration.  The  optimum  design,  shown  on  Table  7,  is 
obtained  a|  the  23rd  iteration.  The  peak  acceleration  is 
342  in/sec^  and  the  cost  is  1.08  watts.  However,  the  comfort 
and  safety  on  each  road  condition  are  much  improved.  As  a 
measure  of  comfort,  the  history  of  absorbed  power  constraint  on 
rough  highway  is  shown  in  Fig.  14.  Cost  and  maximum  violation 
are  plotted  in  Fig.  15. 

The  general  tendency  in  design  optimization  is 

i)  to  increase  stiffness  and  damping  in  the  front  wheel 
suspension  and 

ii)  to  reduce  stiffness  and  damping  in  the  rear  wheel  suspen¬ 
sion. 

This  indicates  that  more  energy  should  be  dissipated  by  the 
front  wheel  suspension  system  to  satisfy  the  given  constraints. 


Table  7.  Initial  and  Optimum  Designs 


Design  Variable  Description 

Initial 

Design 

Optimum 

Design 

Driver  seat  spring  constant  [lb/in] 

b1 

100. 

126.258392 

Spring  constant  of  front 
suspension 

[lb/in] 

b2 

300. 

356.348755 

Spring  constant  of  rear 
suspension 

[lb/in] 

b3 

300. 

274.080944 

Driver  seat  damping 
coefficient 

[lb- sec/in] 

b4 

10. 

2.466310 

Damping  coefficient  of 
front  suspension 

[lb-sec/in] 

b5 

25. 

49.873016 

Damping  coefficient  of 
rear  suspension 

[lb-sec/in] 

b6 

25. 

14.904079 
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ION/  COST 


CONCLUSION 


The  final  results  of  the  optimization,  shown  in  Fig.  15  show 
that  the  absorbed  power  over  the  smooth  highway  has  been 
increased  slightly,  but  that  potential  hazardous  characteris¬ 
tics  of  the  initial  design  have  been  brought  under  control  by 
virtually  eliminating  violations  in  constraints. 

The  success  of  the  method  for  this  simple  dynamic  system  demon¬ 
strates  the  potential  that  exists  in  optimizing  design  for  more 
complex  and  realistic  systems.  The  computer  optimization 
method,  using  the  adjoint  variable  technique  has  been  demon¬ 
strated  to  be  a  useful  tool  in  the  design  of  vehicle  suspen¬ 
sion. 
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ABSTRACT 

A  method  is  presented  for  transient  dynamic  analysis  of  mechan-^ 
ical  systems  composed  of  interconnected  rigid  and  flexible 
bodies  that  undergo  large  angular  displacements.  Gross  dis¬ 
placement  of  elastic  bodies  is  represented  by  superposition 
of  local  linear  elastic  deformation  on  nonlinear  displacement 
of  body  reference  coordinate  systems.  Flexible  bodies  are  thus 
represented  by  combined  sets  of  reference  and  local  elastic 
variables.  Modal  analysis  and  substructuring  of  the  local 
elastic  system  is  employed  to  identify  all  modes  and  eliminate 
the  insignificant  ones.  Equations  of  motion  and  constraints 
are  formulated  in  terms  of  a  minimal  set  of  modal  and  reference 
generalized  coordinates,  which  are  then  dynamically  adjusted  to 
a  time  and  inertia-variant  eigenspectrum. 

Control  system  differential  equations  are  formulated  and 
coupled  with  mechanical  system  models.  The  composite  nonlinear 
system  of  equations  is  then  numerically  integrated  in  the  time 
domain,  to  obtain  complete  system  dynamic  response.  Example 
problems  are  presented  to  demonstrate  the  algorithms  and  com¬ 
puter  code . 

1 .  INTRODUCTION 

Dynamic  analysis  of  large  scale  spatial  nonlinear  elastic 
mechanical  and  control  systems  require  efficient  numerical 
methods.  Modal  techniques  are  desirable  because  they  often 
allow  reduction  of  problems  with  thousands  of  degrees  of  free¬ 
dom  to  manageable  sizes.  The  main  difficulty  with  analyzing 
interconnected  flexible  mechanical  systems,  is  that  composite 
systems  are  generally  inertia-variant  and  their  corresponding 
eigenspectrums  are  time  dependent.  This  necessitates  periodic 
re-evaluation  of  eigenvectors,  an  expensive  process  that  should 
be  avoided. 

To  circumvent  this  problem,  elastic  properties  of  bodies  are 
characterized  locally,  relative  to  body-fixed  reference  frames. 
When  elastic  displacements  are  small,  which  is  often  the  case, 
they  can  be  efficiently  represented  locally  by  modal  analysis 
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techniques  because  the  local  eigenspectrums  are  essentially 
constant.  Thus  it  is  possible  to  efficiently  reduce  the 
dimension  of  a  problem  locally  and  then  transform  the  much 
smaller  problem  to  the  global  system  and  solve  it  using 
existing  numerical  methods. 

The  study  of  system  flexibility  has  been  a  subject  of  major 
interest  to  many  researchers.  In  the  field  of  mechanics ,  there 
have  been  a  number  of  attempts  to  arrive  at  general  algorithms 
for  dynamic  analysis  of  mechanisms  with  elastic  links  [1-11 ]• 
However,  most  of  these  techniques  assume  that  elastic  deforma¬ 
tions  do  not  significantly  affect  gross  body  motion.  Inertia 
forces  of  gross  rigid  body  motion  are  first  determined  and  then 
introduced  as  externally  applied  forces  to  the  elastic  model. 

In  addition,  most  of  these  techniques  are  suitable  for  only 
certain  classes  of  problems. 

Some  investigators  [12-13]  have  solved  simultaneously  for 
coupled  gross  body  and  elastic  deformation.  However,  because 
of  the  highly  nonlinear  nature  of  these  equations,  these  tech¬ 
niques  have  not  utilized  coordinate  reduction  techniques,  nor 
are  they  suitable  for  large  scale  systems. 

A  general  method  for  dynamic  analysis  of  large  scale  inertia- 
variant  systems  has  been  developed  [14-15]  which  utilizes  co¬ 
ordinate  reduction  techniques  commonly  employed  in  structural 
dynamics  and  solves  simultaneously  for  gross  body  motion  and 
elastic  deformation.  The  method  in  [14]  is  applied  to  large 
scale  control  systems  that  contain  elastic  components*  Mechan¬ 
ical  systems  are  considered  as  collections  of  subsystems  called 
bodies,  substructures  or  components.  Finite  element  or  experi¬ 
mental  methods  can  be  used  to  characterize  the  elastic 
properties  of  each  deformable  body.  Energy  equations  of  sub¬ 
elements  are  written  separately  and  the  elements  of  each  body 
(substructure)  are  then  assembled  using  a  Boolean  matrix 
approach.  The  degrees  of  freedom  of  each  body  are  then  reduced 
using  a  component  mode  technique.  Adjacent  substructures  are 
then  interconnected  using  a  Lagrange  multiplier  technique. 

Then  a  coordinate  partitioning  method  [16]  is  employed  to 
eliminate  excess  dependent  equations  of  motion  resulting  from 
imposition  of  the  constraint  equations. 

2.  GENERALIZED  COORDINATES  AND  ENERGY  EQUATIONS 

In  order  to  specify  the  configuration  of  a  body  or  substructure 
it  is  necessary  to  define  a  set  of  generalized  coordinates  such 
that  the  global  position  and  orientation  of  every  infinitesimal 
volume  on  the  body  is  determined  in  terms  of  these  generalized 
coordinates.  As  shown  in  Fig.  1  let  the  XYZ  Cartesian  coordi¬ 
nate  system  represent  an  inertial  frame  and  the  X^’-Z1  axes 
represent  a  Cartesian  coordinate  system  rigidly  attached  to 
some  infinitesimal  volume  on  the  ith  body.  Using  finite 
elements,,  the  ith  body  is  divided  into  a  number  of  rigidly 
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interconnected  subelements .  The  location  of  an  arbitrary 
infinitesimal  volume  on  the  jth  element  of  this  body  is  defined 
in  terms  of  two  sets  of  generalized  coordinates.  The  first  set 
are  reference  coordinates  locating  the  position  of  the  body- 
fixed  coordinate  system  relative  to  the  inertial  XYZ  frame. 

The  second  set  are  elastic  coordinates  characterizing  elastic 
deformation  of  the  body.  Elastic  coordinates  represent  rela¬ 
tive  translational  and  angular  displacement  of  infinitesimal 
volumes  at  nodal  points  on  the  body.  In  addition  the  location 
and  angular  orientation  of  every  infinitesimal  volume  in  each 
element  can  be  approximated  in  terms  of  its  elastic  coordinates 
and  its  shape  function  [17]. 


Let  X1JY1JZ1J>  as  shown  in  Fig.  1 ,  be  a  Cartesian  coordinate 
system  with  its  origin  affixed  to  an  infinitesimal  volume  at 
some  point  on  the  jth  element  of  the  ith  body,  that  rotates 
with  this  infinitesimal  volume.  The  location  of. an  arbitrary 
infinitesimal  volume,  identified  by  the  point  pij  on  this  .  .  . 
element,  can  be  determined  by  specifying  the  position,  of  X^-Y^Z3- 
and  the  location  of  the  point  p^-J  with  respect  to  X1Y1Z1.  Let 
R1  and  and  61  be,  respectively,  the  translational  an<j  .rotation¬ 
al  coordinates  of  X^^-Z1  with  respect  to  XYZ... Let  X'S^Y’-JZ1-!  be 
a  coordinate  system  parallel  to  the  element  X^-JY^Z-'-J  cpQrdi- 
nate  system  and  located  at  the  origin  of  X^Y^Z3-.  Let  e1J  be 
the  nodal  coordinates  of  the  ijth  element  with  respect  to 
X3-!  Y*-J  .  Let .  u^J  .  denote  the  location  of  point  p3-^  with 
respect  to  X^JY^Z^J  .  The  vector  u^J  can  be  written  as 


(D 


where  <(> XJ  is  the  shape  function 
of  the  ijth  element..  The 
position  of  point  pxJ  is  specif¬ 
ied  by 


r^  =  R1  +  A1  Cij  Cijeij 

-p  ~ 

(2) 

where  A1(91)  is  the  transforma¬ 
tion  matrix  from  X^Y’-Z1  to.  XYZ, 
e3-!  is  the  vector  of  elastic 
generalized  coordinates  of  the 
ith. body  defined  relative  to 
Xiyiz3-,  and  C3-^  and  C3-!  are 
transformation  matrices  from 
XijYijz3-!  to  X^Z3-. 


Figure  1 .  Generalized 
Coordinates 


x 


If  the  rotation  of  X1JY3-JZ1J  with  respect  to  X1Y3-Z1  is  small, 
which  is_the  case  for  small  elastic  deformation,  the  matrices 
C3-!  and  C3-J  are  approximately  constant  and  Eq.  2  becomes 
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(3) 


-  R1  +  A1  e1^ 

“P 

where  N^Cx^j,  yij  ,  z*-J)  is  the  modified  shape  function  of  the 
ijth  element  given  by 

Nlj  .  ciJ  gij  (4) 

The  kinetic  energy  of  the  ijth  element  is  obtained  from, the 
velocity  vector  of  the  infinitesimal  volume  at  point  pij  on  the 
element  as 


R  +  A 


+  A  N 


1J 


(5) 


where  (•)  denotes  differentiation  with  respect  to  time.  The 
second  term  of  Eq.  5  can  be  expressed  as 


A1  N« 


in  order  to  isolate  the  velocity  terms.  Equation  5  is  then 
expressed  as 


kV3 


.T 


[r1 


.T 


e1J  ]J 


(7) 


where  I  is  an  identity  matrix.  The  kinetic  energy  of, the  ijth 
element  is  obtained  by  integrating  over  the  volume  V^-J 


/  dV1-) 

il  ~P  ~p 

v1J 


(8) 


.  *  *  •  • 

where  the  mass  density  of  the  ijth  element  at  p1J  is  pxJ  and 
the  superscript  T  implies  transpose  of  a  matrix.  Using  vector 

notation  q1-'  =  [R1  J1  e1-*  ]  and  substituting  Eq.  7  into 
Eq.  8,  kinetic  energy  becomes 


.T 


:ij  ,  Jl 


(9) 
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where  M1J ,  the  mass  matrix  of  the  ijth  element,  is 


M1J  =  /  p 

vij 


iJ 


I 

.  .T 
B1-* 

.  .T  .T 
A1 


.  .  T  .  . 


B1J 

.  .  T  .  . 
B1J  B1J 

T  T 
A1  B^ 


aVJ 


b±JTa1w;LJ 

Nlj  N1^ 


dV 


1J 


(10) 


The  submatrix  N1-'  NX~*  associated  with  element  elastic  general¬ 
ized  coordinates  is  constant  (see  Eqs.  1  and  4). 

Kinetic  energy  T1  of  the  ith  body  is  obtained  by  summing  ov$r 
each  of  its  elements.  Denoting  the  number  of  elements  by  n1 


T1  =  l 

j-1 

-i  ii T  i2T  iniT  T 

In  vector  notation  £  =  [£  ,  £  ,...,£  ]  , 

M1  =  diag  (M11 ,  Ml2  ,  Min  ),  and  Eq.  11  becomes 


(ID 


1  . • T  ... 
T1  =  i  £  M  £ 


(12) 


_i 


The  vector  q  is  partitioned  as 


•  .  •  .T  •  .T  T 

11  -  tii  .  ii  ] 

.T  .T  T 

where  q1  =  [R1  ,  01  ]  are  reference  coordinates  locating 

rp  rp  J  rp  rp 

•  •  •  •  •  4  JL  •  S)  J. 

X1Y1Z1  relative  to  XYZ  and  q^  =  [e1  ,  e1  ,  . . . ,  e 
elastic  generalized  coordinates  of  the  ith  body. 


(13) 


.  iT  T 

]  are 


Since  q|  is  defined  relative  to  X1Y'LZ1  a  constant  Boolean 

matrix  can  be  employed  to  impose  constraints  between  adjacent 
elements  on  the  ith  body.  These  constraints  in  matrix  form 
are 
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I1  -  b|  a1 


(H) 


Differentiating  Eq.  14  with  respect  to  time  and  substituting 
into  Eq.  12  yields 


mi  _  1  *i  Mi  *i 
T  =  2  £  M  3. 


(15) 


where 


.T  .  . 

M1  =  b|  m1  b| 


(16) 


The  kinetic  energy  of  the  ith  body  can  be  partitioned  as 

T 


m1  (q1) 
rrv-i  ' 

i  .  is 

tnfr(q  ) 


Ul 


mrf^L  ' 


m 


ff 


(17) 


where  subscripts  r  and  f  denote  reference  and  flexible  coordi- 

i 

nates  respectively.  Observe  that  m  is  a  constant  matrix  and 

also  that  coupling  between  elastic  and  reference  coordinates 

exists  through  the  nonlinear  matrices  m  and  m  .  Elastic  and 

fr  rf 

reference  coordinates  can  be  dynamically  decoupled  by  neglect- 

ij 

ing  these  two  matrices  and  the  dependence  of  B^  on  elastic 
coordinates. 

Strain  energy  of  the  ijth  elastic  element  is  given  by  [17] 

.  .T 


U 


ij 


7  J.  .  e 


/  eXJ  a dV ^ 


(18) 


where  and  are  the  respective  strain  and  stress  com¬ 

ponents  at  the  infinitesimal  volume.  The  stress-strain 
equation  is  given  by  the  elastic  constitutive  relation 


108 


GX* 


(19) 


* 


t 


and  the  strain-displacement  relation  can  be  expressed  as 
eij  =  Dijuij 


=  D1JN1Je1J  (20) 

where  DLJ  is  a  differential  operator  relating  strains  and  dis¬ 
placements.  Substituting  Eqs.  19  and  20  into  Eq.  18  yields 


U 


ij=7  /  eij 


2  £ 


ijVje^ 


(21) 


where  is  the  ijth  element  stiffness  matrix.  The  element 

stiffness  matrix  is  constant  and  the  strain  energy  expression 
is  independent  of  reference  coordinates  because  elastic, 
generalized  coordinates  are  defined  with  respect  to  X^Y^Z*-. 
The  total  strain  energy  of  the  ith  body  is 


U 


i 


where  =  diag 

Substituting  Eq.  14 


.T  . 

-1  si  -1 

£f  Kff 


0 

0 

»in^ 

•  •  ,Kkff 


) 


into  Eq.  22  and  partitioning  gives 


(22) 


K 


0 

i 

ff  J 


(23) 
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where  q^  and  q^  are  as  defined  before  and  is  the  stiffness 

r  ^f  ff 

matrix  associated  with  the  ith  body  elastic  coordinates. 
Finally  virtual  work  of  the  ijth  element  is 


6W  =  Sij  «£ij 


(24) 


where  is  the  vector  of  generalized  forces  associated  with 

the  generalized  coordinates  £LJ .  All  forces  except  wovkless 
constraint  forces  between  elements  are  included  in  6Vr- J ,  The 
virtual  work  of  the  ith  body  is 


.T  .T  .T  .T  T 

6W1  =  faj  £  }  [  6£^  ] 


(25) 


where  Q  and  Q  are  respectively  generalized  forces  associated 
with  reference  and  elastic  generalized  coordinates. 


3.  SUBSTRUCTURE  SYSTEM  EQUATIONS  OF  MOTION 

The  composite  vector  of  all  system  generalized  coordinates  is 

1T  2T  NT  T 

designated  as  £  =  [£  ,  £  ,...,£  ]  ,  where  N  is  the  total 

number  of  bodies  (substructures)  in  the  system.  Nonlinear 
constraint  equations  between  adjacent  substructures  can  be 
written  in  vector  form  as 

£(£,  t)  »  0  (26) 

where  $(£,  t)  =  [$-|(£,  t) , . . . , *m(£,  t)  and  all  equations  are 
assumed  to  be  independent.  The  equations  of  motion  of  the  ith 
substructure  can  be  written  as  [18]. 


(T*i)T  ‘  (Tli)T  +  fU1.)1  "  S1  +  £Ti 2  -  0  (27) 
£  £  £  £ 


where  \  is  the  vector  of  Lagrange  multipliers. 

Recall  that  all  elastic  deformation  in  the  ith  body  is  defined 
relative  to  the  body  reference  frame  attached  to  node  £  which 
implies  no  elastic  displacement  at  node  £.  Thus 


(28) 


£f  -  0 
£ 


where 


is  the  vector  of  elastic 


.T  .T  T 

node  £.  Denoting  £^  ]  , 

£ 

elastic  generalized  coordinates  of 
£,  Eq.  27  reduces  to 


generalized  coordinates  at 

where  £^  is  the  set  of 
the  ith  body  excluding  node 


1 


M1  (i1)  £X  +  K1^1  =  £, 


t)  +  - 


(29) 


which  is  the  general  form  of  the  equations  of  motion  for  the 
ith  body.  The  vector  F1  absorbs  quadratic  velocity  terms  of 
Eq.  27.  These  equations  along  with  Eq.  26  form  the  constrained 
equations  of  motion  of  the  substructure.  A  coordinate  reduc¬ 
tion  technique  is  now  employed  to  eliminate  insignificant 
elastic  degrees  of  freedom. 


4.  COORDINATE  REDUCTION 

Efficient  solution  of  the  system  equations  of  motion  requires  a 
transformation  from  the  space  of  system  nodal  elastic  general¬ 
ized  coordinates  to  the  space  of  system  modal  generalized  co¬ 
ordinates  with  lower  dimension.  This  transformation  is  a  con¬ 
stant  mapping  [19].  The  problem  is  complicated  by  a  time 
variant  system  transfer  function  which  implies  that  monitoring 
the  frequency  content  of  external  forcing  functions  alone  is 
usually  not  sufficient  for  predicting  mode  excitation  because 
the  mass  matrix  depends  on  generalized  coordinates.  Accord¬ 
ingly,  no  judgement  can  be  made  beforehand  of  the  number  of 
significant  modes  to  be  retained  for  an  accurate  solution. 

The  method  developed  here  is  based  on  solving  the  eigenvalue 
problem  of  the  substructure  once.  From  Fourier  analysis  of  the 
forcing  functions,  an  initial  estimate  of  the  number  of  modes 
to  be  retained  is  made,  and  during  the  simulation  additional 
eigenvectors  are  recalled  or  deleted  as  required.  For  the 
purpose  of  determining  eigenvalues  and  eigenvectors,  if  a  sub¬ 
structure  is  assumed  to  vibrate  freely  about  a  reference  con¬ 
figuration,  Eq.  29  yields 


-1  -1  ,  rl  -1  A 

mff£f  +  -  0 


(30) 
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where  m 
matrices 


and  k^  are  the  respective  mass  and  stiffness 
ff 

associated  with  the  elastic  coordinates.  The  stiff 


ness  matrix  k  is  positive  definite  because  the  reference  co¬ 
ordinate  system  is  fixed.  Equation  30  yields  a  set  of  eigen¬ 
vectors  and  a  modal  matrix.  A  coordinate  transformation  from 
the  physical  elastic  coordinates  to  the  modal  coordinates  is 
obtained  by 


b;  i 


(31) 


where  B*  is  the  modal  matrix  consisting  of  the  eigenvectors 
2 

obtained  from  Eq.  30,  and  y}  is  a  vector  containing  the  modal 
coordinates.  Using  Eq.  31  the  reference  and  elastic  general¬ 
ized  coordinates  are  written  in  terms  of  the  reference  and 
modal  coordinates  as  follows 


I 

0 


(32) 


Substituting  Eq.  32  into  Eq.  29  and  premultiplying  by  the 
transpose  of  the  matrix  in  Eq.  32  yields 


M1 (£1)£1  =  PX(E  >£  >  t)  -  iTi  (2.  t)!  <33> 

2 

The  constraint  Jacobian  matrix  of  Eq.  33  is  evaluated  in  terms 
of  physical  coordinates  and  converted  to  modal  coordinates 
using  Eq.  32  as 

•  • 

!ni=  i_i(9£/8£) 


(34) 


"  ±-i 


I 

L0 
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Equation  33  represents  the  reduced  system  of  equations  for  the 
ith  substructure.  The  number  of  equations  in  Eq.  33  depends  on 
the  number  of  elastic  degrees  of  freedom  required  to  achieve 
the  desired  accuracy. 

5.  CHARACTERIZATION  OF  PLANAR  ELASTIC  SYSTEMS 


Figure  2  illustrates  a  typical  jth  element  denoted  as  ijth  on. 
the  ith  body.  The  element  xfj  axis  is  located  at  an  angle 
relative  to  the  body- fixed  X-*-  coordinate  axis.  The  reference 
coordinates  for  this  body  are  the  translational  variables 
(x^,  yi)  and  the  rotational  variable  0-*-.  This  set  of  coordi¬ 
nates  defines  the  position  of  the  body- fixed  coordinate  system 
X^yi  relative  to  the  inertial  X-Y  frame.  The  elastic  general¬ 
ized  coordinates  are  defined  initially  with  respect  to  a  co¬ 
ordinate  Xij  Y^J  system  which  is  always  parallel  to  the  element 
Xij  Y^-J  axis  with  origin  located  at  the  point  01.  This  set  of 


generalized  coordinates  is  denoted  by  e  ^  (k  =  1 , 


6). 


These  coordinates  are  called  the  nodal  coordinates  and  repre¬ 
sent  the  location  of  the. nodes  and  slopes  of  reference  lines  at 
the  nodes  relative  to  X^j  Y^-j  . 


The  location  of  ap . arbitrary  point  p-k)  on  the  ijth  element, 
with  respect  to  X^-J  Y1^  ,  can  be  expressed  as  [20]: 


(35) 


where 


1  -  K 


3-J 


0 


0 


Lo  1  -  3Ulj  )2+2(?lj  )3  [5lj-2(5lj  )2+(Slj  )3  ] 


•  tj 


0 


0 


0  3(clj  )2-2(clj  )3  )3-  (S1^  )2]  J  (36) 


and 


^  =  x^/^J 


(37) 
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A  transformation  is  employed  to 
define  with  respect  to  X*-  Y* 
as  follows 


w^  *  Cij  e1J 


(38) 


where 

“cos 

6lj 

-sin 

BiJ- 

cij  - 

_sin 

cos 

Blj. 

(39) 


Figure  2.  Two  Dimensional 
Beam  Element 


The  compatibility  conditions  between  elements  on  a  given  sub¬ 
structure  are  simpler  if  the  generalized  elastic  coordinates 
are  defined  with  respect  to  the  body-fixed  coordinate  system  of 
the  substructure.  This  is  accomplished  by  the  following  trans¬ 
formation 


e1J  _  gij  eU 


(40) 


where  e1^  is  the  set. of  generalized  elastic  coordinates  defined 
with  respect  to  X1-  Y1,  and 


C1J 


c1J 

C1 


Lo 


cij 

U1  ■ 


(41) 


with 


tf-3 


cos  3^ 

«  • 

sin  31^ 

o" 

aij 

-sm  3  J 

cos  3^ 

0 

(42) 

0 

0 

1_ 

40  into  Eq. 

38  yields 
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(43) 


w1J  =  C l-J  ^  C1J  eXJ 


It  is  important  to  note  that  the  variatiqn  in  the. 
assumed  to  be  small,  thus  the  matrices  C1!  and  C-'-J 
constant.  Finally,  Eq.  43  can  be  written  as 


w1J  =  N1J 


where 


N1J  =  C1J  C1J 


ij  lj 

The  position  vector  r  of  an  arbitrary  point  p 

P 


element  can  be  expressed  as 


r^  =  R1  +  A1  Nij  e1J 
-P  ~  - 


where 


Di  r  i  iiT 
R  =  fx  y  j 


and 


A1  = 


cos  e 


sin  6' 


-sin  9 1 


cos  0  J 


Differentiating  Eq.  46  with  respect  to  time  gives 


•ii  *i  ,  ;i  „ij  l]  .  Ai  .Tij  »ij 
rJ=R  +  A  N  J  e  J  +  A  NJeJ 
-P 


where 


angle  is 
are  assumed 

(44) 

(45) 

on  the  ijth 

(46) 


(47) 

(48) 
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•i  *i 

A  =  6 


‘-sin  6 


cos  6' 


•i  i' 
=  e1  A1 


ci  _ 
-cos  e 

-sin  6^ 


(49) 


and  A^'  is  the  partial  derivative  of  A*-  vith  respect  to  the 
reference  rotational  degree  of  freedom  01.  Substituting  Eq.  49 

•  ij 

into  Eq.  48  and  writing  r  in  partitioned  form  yields 

”P 

jf  \ 

&  =  [I  A1^1^  aV^]  \  81  (  (50) 

\  e  J) 

Equation  50  is  equivalent  to  Eq.  7  in  the  general  formulation 
with  the  matrix  B-^-j  given  by 


B 


-  A1  N1^ 


(51) 


The  kinetic  energy  expression  for  the  ijth  element  is  given  by 


,lj  -  i  /  Pij  rij  r^  dVij 


V 


ij 


1  *ijT 

=  7  a 


(52) 


where  is  the  element  volume,  p1J 


material  c[ 


T  T 

1J  -  [R1  e1  e1J  ]T  and 


is  the  density  of  element 
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M1J  =  / 


B1^ 


B1J 

B1^  Bij 


1J 


A1  N1J 


.  .T  . 
B^  A1 


N1^  NLj 


dV^ 


(53) 


The  mass  matrix  m  associated  with  the  flexible  coordi- 

ff 

nates  is  the  same  matrix  that  occurs  frequently  in  the  finite 
element  approach  [20]. 

Utilizing  orthonormality  of  A*-',  the  central  term  of  Eq.  53 
becomes 


e1^  J  p^N^V^Y1^ 


/  p1^  BXWJ 


-  eijTm^f  eij 


(54) 


Similarly 

/  p^B^dV1^  =  A1'  /  p^N^dV^e1-3  =  A1 '  Clj  Clj  elj 


V 


ij 


V 


3-J 


=  A1  S1^ 


(55) 


where 


m 

12 


iJ 


6  0  0  6  0 


L0  6  *1J  0  6 


tXh 


(56) 


where  m^J  and  are  the  mass  and  length  of  the  ijth  element 
respectively. 
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The  matrix  product  A 


i 


,T 

A 


i 


is  skew  symmetric, 


i.e. 


» 


L-i  oj 

Thus  from  Eq.  51 

/  pX^B^  A^N^dV1^  *  e1^  f  p^N1^  f  N^dV1^ 

Clj  - 


(57) 


(58) 


where  S  is  a  skew  symmetric  matrix  given  by 


S1J  = 
b1 


mLJ 

w~ 


0 

21 

3*ij 

0 

9 

-2*ij 

-21 

0 

0 

-9 

0 

0 

-3*ij 

0 

0 

•r-l 

c* 

CM 

1 

0 

0 

0 

9 

2£iJ 

0 

21 

-3*ij 

-9 

0 

0 

-21 

0 

0 

2nij 

0 

0 

3£X^ 

0 

0 

(59) 


It  is  noteworthy  that  the  flexibility  mass  matrix  m  is  given 

in  the  finite  element  literature  [20],  thus  it  is  necessary . po 
carry  out  only  the  integrations  required  for  the  matrices  S1^ 

and  S1J  to  completely  define  the  mass  matrix  in  Eq.  53  and 
accordingly  the  kinetic  energy  expression  in  Eq.  52.  Using  the 
above  notations,  the  mass  matrix  is  written  as 
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'h 

. .T  . .T  . ,T 

IJ  qlj  A1 

b  A  (1x2) 

eijTm^fe: 

.  .T  .T 
^  A1 

...  .T  .  . 
S1J 

(symmetric) 


(60) 


(6x2) 


(6x1) 


(6  x6)“ 


where  I  is  a  (2x2)  identity  matrix.  The  total  kinetic  energy 
of  the  ith  body  is  given  by 


T1  =  l  T1-^ 

j-1 


T 

1  *i  „i*i 

7  £  M  £ 


(61) 


where  £1  = 


.  T  .  .T 

[R1  61  ]  represents  the  generalized  coordinates 


of  the  ith  substructure,  c[  the  elastic  coordinates, 


(symmetric) 


i  iT  iT  i’T 
m1  -  sl  A1 


.T  . 
“ff^f 


(62) 


.T  .T 
S1  A1 


„.T  . 

c1  r.1 

S  £f 


m^  is  the  total  mass  of  the  substructure 
S1  is  the  composite  matrix  of  the  S1^  matrices  of  the 
individual  elements 

i  ij 

m  is  the  composite  matrix  of  the  m  matrices  of  the 
ff  ff 

individual  elements 

~i  “ij  . 

S  is  the  composite  matrix  of  the  S  matrices  of  the 


individual  elements 
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It  is  important  to  note  that  each  of  the  above  matrices  is  con¬ 
stant  and  the  implied  assembly  need  only  be  done  once.  This  is 
because  the  elastic  generalized  coordinates  of  each  substruc¬ 
ture  are  defined  locally  which  leads  to  time  invariant  com¬ 
patibility  conditions  between  the  elements. 

Neglecting  shear  deformation,  the  strain  energy  of  the  jth 
element  on  the  ith  body  is  given  by 


U 


0 


)2  +  EijAiJ(uiJ,)2]dxij 


(63) 


where  primes  indicate  derivatives  with  respect  to  x1  j  ,  E1^  is 
th$  modulus  of  elasticity,  I1J  is  the  second  area  moment  and 
is  the  element  cross-sectional  area. 


Using  Eq. 


35,  Eq.  63  becomes 

’  *T«a  eij 


cijT  eij  eij 


_  i  Kij 


^ff  ^ 


(64) 


where  is  the  element  stiffness  matrix  defined  with  respect 

ff 

to  the  body-fixed  coordinate  system. 

Again,  the  stiffness  matrix  is  constant  because  the  elastic 
generalized  coordinates  are  defined  with  respect  to  the  body- 
fixed  coordinate  system.  The  total  strain  energy  of  the  ith 
body  is  obtained  as  follows 


n 


ui  -  X  ulJ 


‘0 

0 


0 


KffJ 


1  n1 Vn1 

=  j  a  K  £ 


(65) 
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where  K 


i 

ff 


is  the  assembled  stiffness  matrix  of  the  ith  body  and 


K1  = 


0 

LO 


0 


KifJ 


(66) 


6.  NUMERICAL  EXAMPLE 

The  Peaucillier  mechanism  shown  in  Fig.  3  is  designed  to 
generate  a  straight  line  path.  Its  geometry  is  such  that 
BC  -  BP  =  EC  =  EP  =  0.359m  and  AB  =  AE  =  0.254m.  Points  A,  C 
and  P  should  always  lie  on  a  straight  line  passing  through  A. 
The  mechanism  always  satisfies  the  condition  AC  x  AP  =  C,  where 
C  is  a  constant  called  the  inversion  constant.  In  case 
AD  =  CD,  point  C  must  trace  a  circular  arc  and  point  P  should 
follow  an  exact  straight  line.  However  this  will  not  be  the 
case  when  flexibility  of  the  links  is  considered. 

All  mechanism  members  are  assumed  constructed  of  steel  bars  of 
circular  cross-sectional  area  of  diameter  0.01m.  The  input 
crank  OC  is  assumed  to  rotate  such  that 

0  =  sin  200t 

where  0  is  the  angular  rotation  of  the  input  link  with  respect 
to  the  XY  inertial  frame  and  t  is  time. 

Links  BP  and  EP  are  assumed  flexible  and  their  body- fixed  co¬ 
ordinate  systems  are  located  at  their  corresponding  centers  of 
mass.  Each  of  these  links  are  divided  into  two  beam  elements. 
Because  of  the  flexibility  of  these  links,  point  P  will  no 
longer  move  in  a  straight  line.  The  deviation  $  of  point  P  is 
plotted  in  Fig.  4  versus  £2,  where  £2  =  200 1.  The  figure  shows 
two  solutions.  The  first  is  a  two-mode  solution,  in  which  the 
deformation  of  each  flexible  link  is  described  by  its  corres¬ 
ponding  two  lowest  modes.  Second,  a  four-mode  solution  des¬ 
cribed  by  the  lowest  four  modes  of  each  flexible  link.  The 
figure  indicates  that  the  first  two  modes  of  each  link  are 
dominant.  The  displacement  3  can  be  substantially  reduced  if 
the  amplitude  of  the  dominant  modes  is  kept  small. 

Equation  (33) ,  representing  the  equations  of  motion  of  the  ith 
substructure,  is  expressed  in  terms  of  reference  and  modal  co¬ 
ordinates.  Therefore  a  modal  control  method  [21]  can  be 
employed  to  control  the  deviation  0  of  point  P.  Equation  (33), 
in  this  case,  can  be  written  as 


121 


Figure  3.  Peaucillier  Lipkin  Figure  4.  The  Two-Mode  and 
Mechanism  Four-Mode  Solution 


» 


-  P^E.E.t)  -  £Ti(E.t)2  -  u(£,t) 

£ 

Where  u<£,t)  is  a  control  input.  This  control  input  represents 
a  controller  to  the  system  designed  in  such  a  way  to  reduce  the 
deviation  $  of  point  P.  This  displacement  can  be  measured 
physically  using  a  displacement  meter.  However,  8  can  also  be 
expressed  mathematically  as  a  function  of  the  elastic  modes 
which  are  obtained  in  the  solution  of  Eq.  33.  Controlling 
these  modes  will  reduce  the  deviation  of  P  from  the  desired 
path.  In  the  present  example  u(£,t)  is  designed  as  a  propor¬ 
tional  controller  which  is  expressed  mathematically  as 

u(£,t)  -  G  £ 

where  G  is  a  diagonal  matrix  containing  the  gain  factors 
(Fig.  5).  This  matrix  is  defined  such  that  the  dominant  modes 
are  multiplied  by  equal  gain.  The  modal  coordinates,  for 
different  values  of  gain  are  shown  in  Figs.  6  to  9.  Figure  10 
shows  8  of  the  four-mode  solution  for  different  gains.  It  can 
be  seen  that  increasing  the  gain  leads  to  a  significant 
reduction  in  the  modal  coordinates  and,  accordingly,  to  a 
reduction  in  8.  This  increase  in  the  gain  is  physically  equiv¬ 
alent  to  a  change  in  the  elastic  member  stiffness  which  can  be 
accomplished  by  changing  the  member  dimensions  and  shapes. 
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Figure  5.  The  System 
Representation  with 
Proportional  Controller 


Figure  6.  Effect  of  the  Gain 
on  the  First  Mode 
of  Link  BP 


oo 


CO 


n  rad 


Figure  7.  Effect  of  the  Gain 
on  the  Second  Mode 
of  Link  BP 


Figure  8.  Effect  of  the  Gain 
on  the  First  Mode 
of  Link  EP 


Figure  9.  Effect  of  the  Gain 
on  the  Second  Mode 
of  Link  EP 


Figure  10.  Four-Mode  Solu¬ 
tion  for  Different 
Gain  Factors 


7 •  CONCLUSION 

A  method  for  the  dynamic  analysis  of  large  scale  inertia- 
variant  mechanical  systems  with  flexible  components  is 
presented.  The  resultant  method  is  capable  of  analyzing  com¬ 
plex  mechanical  systems  and  it  utilizes  coordinate  reduction 
techniques.  The  final  form  of  the  system  equations  of  motion 
is  expressed  in  terms  of  the  reference  and  modal  coordinates. 
The  applicability  of  the  method  to  control  systems  is  demon¬ 
strated  through  a  simple  example.  The  modal  control  technique 
is  employed  to  control  the  dominant  modes  of  the  system  com¬ 
ponents  . 
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